src/MEDMEM/MEDMEM_Connectivity.cxx

   1 //  MED MEDMEM : MED files in memory
   2 //
   3 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  21 //
  22 //
  23 //
  24 //  File   : MEDMEM_Connectivity.cxx
  25 //  Module : MED
  26
  27 using namespace std;
  28 #include "MEDMEM_Connectivity.hxx"
  29 #include "MEDMEM_Family.hxx"
  30 #include "MEDMEM_CellModel.hxx"
  31
  32 #include "MEDMEM_SkyLineArray.hxx"
  33 #include "MEDMEM_ModulusArray.hxx"
  34
  35 #include "MEDMEM_STRING.hxx"
  36
  37 /*!
  38    Default Constructor. /n
  39    Default for Entity is MED_CELL and type of Connectivity is MED_NODAL */
  40 //--------------------------------------------------------------//
  41 CONNECTIVITY::CONNECTIVITY(medEntityMesh Entity /* =MED_CELL */) :
  42 //--------------------------------------------------------------//
  43                                 _entity(Entity),
  44                                 _typeConnectivity(MED_NODAL),
  45                                 _numberOfTypes(0),
  46                                 _geometricTypes((medGeometryElement*)NULL),
  47                                 _type((CELLMODEL*)NULL),
  48                                 _entityDimension(0),
  49                                 _count((int*)NULL),
  50                                 _nodal((MEDSKYLINEARRAY*)NULL),
  51                                 _descending((MEDSKYLINEARRAY*)NULL),
  52                                 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
  53                                 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
  54                                 _neighbourhood((MEDSKYLINEARRAY*)NULL),
  55                                 _constituent((CONNECTIVITY*)NULL)
  56 {
  57    MESSAGE("CONNECTIVITY(medEntityMesh Entity=MED_CELL)");
  58 }
  59
  60 /*!
  61    Constructor. /n
  62    Default for Entity is MED_CELL */
  63 //------------------------------------------------------------------------------//
  64 CONNECTIVITY::CONNECTIVITY(int numberOfTypes,medEntityMesh Entity /* =MED_CELL */):
  65 //------------------------------------------------------------------------------//
  66                                 _entity(Entity),
  67                                 _typeConnectivity(MED_NODAL),
  68                                 _numberOfTypes(numberOfTypes),
  69                                 _entityDimension(0),
  70                                 _nodal((MEDSKYLINEARRAY*)NULL),
  71                                 _descending((MEDSKYLINEARRAY*)NULL),
  72                                 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
  73                                 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
  74                                 _neighbourhood((MEDSKYLINEARRAY*)NULL),
  75                                 _constituent((CONNECTIVITY*)NULL)
  76 {
  77   MESSAGE("CONNECTIVITY(int numberOfTypes,medEntityMesh Entity=MED_CELL)");
  78   _geometricTypes = new medGeometryElement[numberOfTypes];
  79   _type = new CELLMODEL[numberOfTypes];
  80   _count = new int[numberOfTypes+1];
  81 }
  82
  83 /*!
  84   Copy Constructor.
  85 */
  86 //------------------------------------------------------//
  87 CONNECTIVITY::CONNECTIVITY (const CONNECTIVITY & m):
  88 //----------------------------------------------------//
  89                                 _entity           (m._entity),
  90                                 _typeConnectivity (m._typeConnectivity),
  91                                 _numberOfTypes    (m._numberOfTypes),
  92                                 _entityDimension  (m._entityDimension),
  93                                 _numberOfNodes    (m._numberOfNodes)
  94 {
  95  if (m._geometricTypes != NULL)
  96  {
  97     _geometricTypes = new medGeometryElement[_numberOfTypes];
  98     memcpy(_geometricTypes,m._geometricTypes,_numberOfTypes*sizeof(medGeometryElement));
  99  }
 100  else
 101     _geometricTypes = (medGeometryElement *) NULL;
 102
 103  if (m._type != NULL)
 104  {
 105     _type = new CELLMODEL[_numberOfTypes];
 106     for (int i=0;i<_numberOfTypes;i++)
 107         _type[i] = CELLMODEL(m._type[i]);
 108  }
 109  else
 110     _type = (CELLMODEL *) NULL;
 111
 112  if (m._count != NULL)
 113  {
 114       _count = new med_int[_numberOfTypes+1];
 115       memcpy(_count,m._count,(_numberOfTypes+1)*sizeof(med_int));
 116  }
 117  else
 118     _count = (med_int *) NULL;
 119
 120  if (m._nodal != NULL)
 121     _nodal = new MEDSKYLINEARRAY(* m._nodal);
 122  else
 123     _nodal = (MEDSKYLINEARRAY *) NULL;
 124
 125  if (m._descending != NULL)
 126     _descending = new MEDSKYLINEARRAY(* m._descending);
 127  else
 128     _descending = (MEDSKYLINEARRAY *) NULL;
 129
 130  if (m._reverseNodalConnectivity != NULL)
 131     _reverseNodalConnectivity = new MEDSKYLINEARRAY(* m._reverseNodalConnectivity);
 132  else
 133     _reverseNodalConnectivity = (MEDSKYLINEARRAY *) NULL;
 134
 135  if (m._reverseDescendingConnectivity != NULL)
 136     _reverseDescendingConnectivity = new MEDSKYLINEARRAY(* m._reverseDescendingConnectivity);
 137  else
 138     _reverseDescendingConnectivity = (MEDSKYLINEARRAY *) NULL;
 139
 140  if (m._neighbourhood != NULL)
 141     _neighbourhood = new MEDSKYLINEARRAY(* m._neighbourhood);
 142  else
 143     _neighbourhood = (MEDSKYLINEARRAY *) NULL;
 144
 145  if (m._constituent != NULL)
 146     _constituent = new CONNECTIVITY(* m._constituent);
 147  else
 148     _constituent = (CONNECTIVITY *) NULL;
 149 }
 150
 151 /*!
 152    Destructor./n
 153    desallocates existing pointers */
 154 //----------------------------//
 155 CONNECTIVITY::~CONNECTIVITY()
 156 //----------------------------//
 157 {
 158   MESSAGE("Destructeur de CONNECTIVITY()");
 159
 160   if (_geometricTypes != NULL)
 161      delete [] _geometricTypes;
 162   if (_type != NULL)
 163      delete [] _type;
 164   if (_count != NULL)
 165      delete[] _count;
 166   if (_nodal != NULL)
 167      delete _nodal;
 168   if (_descending != NULL)
 169      delete _descending;
 170   if (_reverseNodalConnectivity != NULL)
 171      delete _reverseNodalConnectivity;
 172   if (_reverseDescendingConnectivity != NULL)
 173      delete _reverseDescendingConnectivity;
 174   if (_constituent != NULL)
 175      delete _constituent;
 176 }
 177
 178 /*!
 179    set _constituent to Constituent
 180    be aware desallocation of _constituent is done by CONNECTIVITY:~CONNECTIVITY
 181    throws an exception if Constituent = MED_CELL
 182    A DOCUMENTER
 183     */
 184 //----------------------------------------------------------//
 185 void CONNECTIVITY::setConstituent(CONNECTIVITY * Constituent)
 186                     throw (MEDEXCEPTION)
 187 //----------------------------------------------------------//
 188 {
 189   medEntityMesh Entity = Constituent->getEntity();
 190   if (Entity == MED_CELL)
 191     throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : could not set constituent on MED_CELL !"));
 192
 193   if ((Entity == MED_EDGE)&(_entityDimension == 3))
 194   {
 195     if (_constituent == NULL)
 196       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : Entity not found !"));
 197     _constituent->setConstituent(Constituent);
 198   }
 199   else
 200     _constituent = Constituent;
 201 }
 202
 203 /*! Duplicated Types array in CONNECTIVITY object. */
 204 //--------------------------------------------------------------------//
 205 void CONNECTIVITY::setGeometricTypes(const medGeometryElement * Types,
 206                                      const medEntityMesh Entity)
 207                                      throw (MEDEXCEPTION)
 208 //--------------------------------------------------------------------//
 209 {
 210   if (Entity == _entity)
 211     for (int i=0; i<_numberOfTypes; i++)
 212     {
 213       _geometricTypes[i] = Types[i];
 214       _type[i] = CELLMODEL(_geometricTypes[i]);
 215     }
 216   else
 217     {
 218     if (_constituent == NULL)
 219         throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setGeometricTypes : Entity not found !"));
 220     _constituent->setGeometricTypes(Types,Entity);
 221   }
 222 }
 223
 224 /*! A DOCUMENTER */
 225 //--------------------------------------------------------------------//
 226 void CONNECTIVITY::setCount(const int * Count, const medEntityMesh Entity)
 227                             throw (MEDEXCEPTION)
 228 //--------------------------------------------------------------------//
 229 {
 230   if (Entity == _entity)
 231   {
 232     int * index = new int[Count[_numberOfTypes]];
 233     index[0]=1;
 234     _count[0]=1;
 235     for (int i=0; i<_numberOfTypes; i++) {
 236       _count[i+1] = Count[i+1];
 237       int NumberOfNodesPerElement = _type[i].getNumberOfNodes() ;
 238       for (int j=_count[i]; j<_count[i+1]; j++)
 239         index[j] = index[j-1]+NumberOfNodesPerElement;
 240     }
 241     // allocate _nodal
 242     if (_nodal != NULL) delete _nodal;
 243     _nodal = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,index[_count[_numberOfTypes]-1]-1);
 244     _nodal->setIndex(index);
 245     delete[] index;
 246   }
 247   else
 248   {
 249     if (_constituent == NULL)
 250       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setCount : Entity not found !"));
 251     _constituent->setCount(Count,Entity);
 252   }
 253 }
 254
 255 //--------------------------------------------------------------------//
 256 void CONNECTIVITY::setNodal(const int * Connectivity,
 257                             const medEntityMesh Entity,
 258                             const medGeometryElement Type)
 259                             throw (MEDEXCEPTION)
 260 //--------------------------------------------------------------------//
 261 {
 262   if (_entity == Entity)
 263   {
 264     // find geometric type
 265     bool find = false;
 266     for (int i=0; i<_numberOfTypes; i++)
 267       if (_geometricTypes[i] == Type) {
 268         find = true;
 269         int NumberOfNodePerElements = _type[i].getNumberOfNodes() ;
 270         //_nodal->setI(i+1,Connectivity);
 271         for( int j=_count[i];j<_count[i+1]; j++)
 272           _nodal->setI(j,Connectivity+(j-_count[i])*NumberOfNodePerElements);
 273       }
 274     if (!find)
 275       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : geometric type not found"));
 276   } else
 277   {
 278     if (_constituent == NULL)
 279       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : Entity not found !"));
 280       _constituent->setNodal(Connectivity,Entity,Type);
 281   }
 282 }
 283
 284 /*! A DOCUMENTER */
 285 //------------------------------------------------------------------------------------------//
 286 void CONNECTIVITY::calculateConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity)
 287 //------------------------------------------------------------------------------------------//
 288 {
 289   MESSAGE("CONNECTIVITY::calculateConnectivity");
 290
 291   // a temporary limitation due to calculteDescendingConnectivity function !!!
 292   if ((_entityDimension==3) & (Entity==MED_EDGE))
 293     throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build EDGE !");
 294
 295   if (Entity==_entity)
 296     if (ConnectivityType==MED_NODAL)
 297       calculateNodalConnectivity();
 298     else
 299       if (Entity==MED_CELL)
 300         calculateDescendingConnectivity();
 301       else
 302         throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build DESCENDING !");
 303   if (Entity!=_entity) {
 304     calculateDescendingConnectivity();
 305     _constituent->calculateConnectivity(ConnectivityType,Entity);
 306   }
 307 }
 308
 309 /*!  Give, in full or no interlace mode (for nodal connectivity),
 310      descending or nodal connectivity.
 311
 312       You must give a <medEntityMesh> (ie:MED_EDGE)
 313       and a <medGeometryElement> (ie:MED_SEG3).  */
 314
 315 //------------------------------------------------------------//
 316 void CONNECTIVITY::updateFamily(vector<FAMILY*> myFamilies)
 317 //------------------------------------------------------------//
 318 {
 319   const char * LOC = "CONNECTIVITY::updateFamily(vector<FAMILY*>) ";
 320   BEGIN_OF(LOC);
 321
 322   int numberOfFamilies = myFamilies.size();
 323   if (numberOfFamilies == 0 ) {
 324     MESSAGE(LOC<<"No family");
 325     return;
 326   }
 327   // does we do an update ?
 328   if ((_constituent != NULL)&(_descending != NULL)) {
 329     MESSAGE(LOC<<"Constituent is already defined");
 330     return;
 331   }
 332
 333   if ((_constituent != NULL)&(_descending == NULL)) {
 334     if (myFamilies[0]->getEntity() != _constituent->getEntity()) {
 335       MESSAGE(LOC<<"Family and constituent entity are different. We do nothing");
 336       return;
 337     }
 338
 339     // well we could go !
 340     CONNECTIVITY * oldConstituent = _constituent;
 341
 342 //      for(int i=0; i<numberOfFamilies; i++) {
 343 //        FAMILY * myFamily = myFamilies[i];
 344 //        MESSAGE(LOC<<"updating the family (BEGIN) : " << *myFamily);
 345 //      }
 346
 347     _constituent = (CONNECTIVITY *)NULL;
 348     // for instance we must have nodal connectivity in constituent :
 349     if (oldConstituent->_nodal == NULL)
 350       {
 351         MESSAGE(LOC<<"We have no nodal connectivity of sub cell");
 352         throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
 353       }
 354     int oldNumberOfFace = oldConstituent->_nodal->getNumberOf();
 355     const int * oldConstituentValue = oldConstituent->_nodal->getValue();
 356     const int * oldConstituentIndex = oldConstituent->_nodal->getIndex();
 357
 358     calculateDescendingConnectivity();
 359
 360     int newNumberOfFace = _constituent->_nodal->getNumberOf();
 361     const int * newConstituentValue = _constituent->_nodal->getValue();
 362     const int * newConstituentIndex = _constituent->_nodal->getIndex();
 363
 364     const int * newReverseDescendingIndex =
 365       _reverseDescendingConnectivity->getIndex();
 366
 367     const int * newDescendingIndex = _descending->getIndex();
 368     const int * newDescendingValue = _descending->getValue();
 369
 370     // loop on all family,
 371     //   for all constituent in family, we get it's old connectivity
 372     //   (with oldCconstituentValue and oldConstituentIndex)
 373     //   and search the constituent in newConstituentValue with class
 374     //   ModulusArry
 375     //
 376     //   when a new face is found, replace old constituent
 377     //   number in family with new one
 378     //   If face is not rigth oriented, we must change _descending attribute
 379     //   and _reverseDescendingConnectivity (see calculateDescendingConnectivity()).
 380
 381     // Voila a toi de jouer Nadir :-)
 382
 383     // First we get the renumbering from the oldCconstituentValue and
 384     // oldConstituentIndex in the the new one, newConstituentValue and
 385     // newConstituentIndex with a negative sign if the face is not
 386     // right orented
 387
 388     int * renumberingFromOldToNew = new int [oldNumberOfFace];
 389     int index1 = 0;
 390     int indexm1 = 0;
 391
 392     _constituent->calculateReverseNodalConnectivity();
 393
 394     for (int iOldFace=0;iOldFace<oldNumberOfFace;iOldFace++)
 395       {
 396         int index = 0;
 397
 398         renumberingFromOldToNew[iOldFace] = iOldFace+1;
 399         //        renumberingFromOldToNew[iOldFace] = 999999;
 400
 401         int face_it_beginOld = oldConstituentIndex[iOldFace];
 402         int face_it_endOld = oldConstituentIndex[iOldFace+1];
 403         int face_size_itOld = face_it_endOld - face_it_beginOld;
 404
 405         const int* NodesLists = oldConstituentValue + (face_it_beginOld-1);
 406         int face_size_itNew;
 407
 408         const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity();
 409         const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex();
 410
 411         // set an array wich contains faces numbers arround first node
 412         int BeginIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]-1];
 413         int EndIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]];
 414         int NumberOfFacesInList=EndIndexFaceArrayFirstNode-BeginIndexFaceArrayFirstNode;
 415
 416         int * FacesList = new int[NumberOfFacesInList];
 417
 418         for (int l=BeginIndexFaceArrayFirstNode; l<EndIndexFaceArrayFirstNode; l++){
 419           FacesList[l-BeginIndexFaceArrayFirstNode]=reverseFaceNodal[l-1];
 420         }
 421         // foreach node in sub cell, we search elements which are in common
 422         // at the end, we must have only one !
 423
 424         for (int nodeFaceOld=1; nodeFaceOld<face_size_itOld; nodeFaceOld++)
 425           {
 426             int NewNumberOfFacesInList = 0;
 427             int * NewFacesList = new int[NumberOfFacesInList];
 428
 429             for (int l1=0; l1<NumberOfFacesInList; l1++) {
 430               for (int l2=reverseFaceNodalIndex[NodesLists[nodeFaceOld]-1]; l2<reverseFaceNodalIndex[NodesLists[nodeFaceOld]]; l2++) {
 431                 if (FacesList[l1]<reverseFaceNodal[l2-1])
 432                   // increasing order : FacesList[l1] are not in elements list of node l
 433                   break;
 434                 if (FacesList[l1]==reverseFaceNodal[l2-1]) {
 435                   // we have found one
 436                   NewFacesList[NewNumberOfFacesInList]=FacesList[l1];
 437                   NewNumberOfFacesInList++;
 438                   break;
 439                 }
 440               }
 441             }
 442             NumberOfFacesInList = NewNumberOfFacesInList;
 443             delete [] FacesList;
 444             FacesList = NewFacesList;
 445           }
 446
 447         if (!NumberOfFacesInList==0)
 448           {
 449             if (NumberOfFacesInList>1)
 450               throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one face found ("<<NumberOfFacesInList<<") ! " <<FacesList[0]<<" "<<FacesList[1] ));
 451
 452             MEDMODULUSARRAY modulusArrayOld(face_size_itOld,NodesLists);
 453
 454             int face_it_beginNew = newConstituentIndex[FacesList[0]-1];
 455             int face_it_endNew = newConstituentIndex[FacesList[0]];
 456             face_size_itNew = face_it_endNew - face_it_beginNew;
 457
 458             const int * newNodesLists = newConstituentValue+newConstituentIndex[FacesList[0]-1]-1;
 459             MEDMODULUSARRAY modulusArrayNew(face_size_itNew,newNodesLists);
 460
 461             int retCompareNewOld = modulusArrayNew.compare(modulusArrayOld);
 462
 463             //SCRUTE(retCompareNewOld);
 464
 465             // Real new face found
 466
 467             if(retCompareNewOld == 1)
 468               {
 469                 renumberingFromOldToNew[iOldFace] = FacesList[0];
 470                 index = 1;
 471                 index1++;
 472               }
 473
 474             // Reverse new face found
 475
 476             if(retCompareNewOld == -1)
 477               {
 478                 renumberingFromOldToNew[iOldFace] = FacesList[0];
 479                 index = 1;
 480                 indexm1++;
 481
 482                 int face_it_begin = newReverseDescendingIndex[FacesList[0]-1];
 483                 int face_it_end = newReverseDescendingIndex[FacesList[0]];
 484                 int face_size_it = face_it_end - face_it_begin;
 485
 486                 if (face_size_it == 1)
 487                   throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
 488
 489                 if (face_size_it > 2)
 490                   throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This face/edge should not be a (d-1) cell because it has "<<face_size_it<<" neighbouring"));
 491
 492                 // we have always 2 neighbourings
 493                 int cell1 = _reverseDescendingConnectivity->getIJ(FacesList[0],1);
 494                 int cell2 = _reverseDescendingConnectivity->getIJ(FacesList[0],2);
 495                 // PROVISOIRE : en attendant que le SKYLINEARRAY de ReverseDescending soit correct (sans le zero)
 496                 //                  if (cell2 == 0)
 497                 //                    throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
 498
 499                 if (cell2 != 0) { // we are not on border !!!!
 500
 501                   _reverseDescendingConnectivity->setIJ(FacesList[0],1,cell2);
 502                   // Updating _constituent->_nodal because of reversity
 503                   const int * oldArray = oldConstituentValue+face_it_beginOld-1;
 504                   for(int iarray=1;iarray<=face_size_itNew;iarray++){
 505                     _constituent->_nodal->setIJ(FacesList[0],iarray,oldArray[iarray-1]);
 506                   }
 507
 508                   // Updating _reverseDescendingConnectivity
 509
 510
 511                   _reverseDescendingConnectivity->setIJ(FacesList[0],2,cell1);
 512
 513                   // Updating _descending for cell1 and cell2
 514                   for(int iface=newDescendingIndex[cell1-1];iface<=newDescendingIndex[cell1];iface++)
 515                     if (_descending->getIndexValue(iface)==FacesList[0])
 516                       _descending->setIndexValue(iface,-FacesList[0]);
 517                     else if (_descending->getIndexValue(iface)==-FacesList[0])
 518                       _descending->setIndexValue(iface,FacesList[0]);
 519                   // else nothing to do
 520
 521                   for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];iface++)
 522                     if (_descending->getIndexValue(iface)==FacesList[0])
 523                       _descending->setIndexValue(iface,-FacesList[0]);
 524                     else if (_descending->getIndexValue(iface)==-FacesList[0])
 525                       _descending->setIndexValue(iface,FacesList[0]);
 526                   // else nothing to do
 527
 528                 } else {// else we are on border and we do nothing !!!!!!!!
 529                   INFOS("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
 530                   INFOS(LOC<<" Boudary FACE "<<iOldFace+1<<" are wrong oriented !");
 531                   INFOS("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
 532                 }
 533               }
 534
 535             if(index == 0)
 536               {
 537                 INFOS(LOC<<"Renumbering problem with the Face connectivity given by the User and the new Connectivity computed");
 538                 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a Face connectivity problem"));
 539               }
 540           }
 541         delete[] FacesList;
 542       }
 543
 544     MESSAGE(LOC<<"The Renumbering is finished and the status is");
 545
 546     // Updating the Family
 547
 548     for(int i=0; i<numberOfFamilies; i++) {
 549       FAMILY * myFamily = myFamilies[i];
 550
 551       MEDSKYLINEARRAY * number = myFamily->getnumber();
 552       int numberOfLines_skyline = number->getNumberOf();
 553       const int * index_skyline = number->getIndex();
 554
 555       for (int i=0;i<numberOfLines_skyline;i++) {
 556         for (int j=index_skyline[i]; j<index_skyline[i+1];j++) {
 557           number->setIndexValue(j,renumberingFromOldToNew[number->getIndexValue(j)-1]);
 558         }
 559       }
 560       MESSAGE(LOC<<"updating the family (END) : " << *myFamily);
 561     }
 562
 563     delete oldConstituent ;
 564     delete [] renumberingFromOldToNew;
 565   }
 566
 567
 568   END_OF(LOC);
 569   return;
 570 }
 571
 572 //------------------------------------------------------------------------------------------------------------------//
 573 const med_int * CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type)
 574 //------------------------------------------------------------------------------------------------------------------//
 575 {
 576   const char * LOC = "CONNECTIVITY::getConnectivity";
 577   BEGIN_OF(LOC);
 578
 579   MEDSKYLINEARRAY * Connectivity;
 580   if (Entity==_entity) {
 581
 582     if (ConnectivityType==MED_NODAL)
 583       {
 584         calculateNodalConnectivity();
 585         Connectivity=_nodal;
 586       }
 587     else
 588       {
 589         calculateDescendingConnectivity();
 590         Connectivity=_descending;
 591       }
 592
 593     if (Connectivity!=NULL)
 594       if (Type==MED_ALL_ELEMENTS)
 595         return Connectivity->getValue();
 596       else {
 597         for (med_int i=0; i<_numberOfTypes; i++)
 598           if (_geometricTypes[i]==Type)
 599             //return Connectivity->getI(i+1);
 600             return Connectivity->getI(_count[i]);
 601         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
 602       }
 603     else
 604       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
 605   } else
 606     if (_constituent != NULL)
 607       return _constituent->getConnectivity(ConnectivityType,Entity,Type);
 608
 609   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 610 }
 611
 612 /*!  Give morse index array to use with
 613      getConnectivity(MED_FULL_INTERLACE,mode,entity,MED_ALL_ELEMENTS).
 614
 615       Each value give start index for corresponding entity in connectivity array.
 616
 617       Example : i-th element, j-th node of it :
 618       - In C mode : Connectivity[ConnectivityIndex[i]-1+j-1]
 619       - In fortran mode : Connectivity[ConnectivityIndex[i]+j] */
 620 //-----------------------------------------------------------------------------------------------//
 621 const med_int * CONNECTIVITY::getConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity)
 622 //------------------------------------------------------------------------------------------------//
 623 {
 624   const char * LOC = "CONNECTIVITY::getConnectivityIndex";
 625   BEGIN_OF(LOC);
 626
 627   MEDSKYLINEARRAY * Connectivity;
 628   if (Entity==_entity) {
 629
 630     if (ConnectivityType==MED_NODAL)
 631       Connectivity=_nodal;
 632     else
 633       Connectivity=_descending;
 634
 635     if (Connectivity!=NULL)
 636       return Connectivity->getIndex();
 637     else
 638       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
 639
 640   } else
 641     if (_constituent != NULL)
 642       return _constituent->getConnectivityIndex(ConnectivityType,Entity);
 643
 644   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 645 }
 646
 647 /*! A DOCUMENTER */
 648 //--------------------------------------------------------------//
 649 const CELLMODEL & CONNECTIVITY::getType(medGeometryElement Type) const
 650 //--------------------------------------------------------------//
 651 {
 652   const char * LOC = "CONNECTIVITY::getType";
 653   BEGIN_OF(LOC);
 654
 655   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
 656     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE !"));
 657   for (med_int i=0; i<_numberOfTypes; i++)
 658     if (_geometricTypes[i]==Type)
 659       return _type[i];
 660   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
 661 }
 662
 663 /*!  Returns the number of elements of type <med geometrie element>.
 664      Note : not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
 665 //------------------------------------------------------------------------//
 666 med_int CONNECTIVITY::getNumberOfNodesInType(medGeometryElement Type) const
 667 //------------------------------------------------------------------------//
 668 {
 669   const char * LOC = "CONNECTIVITY::getNumberOfNodesInType";
 670   BEGIN_OF(LOC);
 671
 672   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
 673     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!"));
 674   for (med_int i=0; i<_numberOfTypes; i++)
 675     if (_geometricTypes[i]==Type)
 676       return _type[i].getNumberOfNodes();
 677   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
 678 }
 679
 680 /*!  Returns the number of geometric sub cells of <med geometrie element> type.
 681 not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
 682 //------------------------------------------------------------------------//
 683 med_int CONNECTIVITY::getNumberOfSubCellInType(medGeometryElement Type) const
 684 //------------------------------------------------------------------------//
 685 {
 686   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
 687     throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!");
 688   for (med_int i=0; i<_numberOfTypes; i++)
 689     if (_geometricTypes[i]==Type)
 690       return _type[i].getNumberOfConstituents(1);
 691   throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement not found !");
 692 }
 693
 694 /*!  Returns the number of elements of type <med geometrie element>.
 695
 696      Note :
 697       - Implemented for MED_ALL_ELEMENTS
 698       - Not implemented for MED_ALL_ENTITIES (A VERIFIER)
 699       - Not implemented for MED_NONE */
 700 //-----------------------------------------------------------------------------------//
 701 med_int CONNECTIVITY::getNumberOf(medEntityMesh Entity, medGeometryElement Type) const
 702 //-----------------------------------------------------------------------------------//
 703 {
 704   const char * LOC = "CONNECTIVITY::getNumberOf";
 705   BEGIN_OF(LOC);
 706
 707   MESSAGE(LOC<<" Entity = "<< Entity << ", _entity = "<<_entity);
 708
 709   if (Entity==_entity) {
 710     if (Type==MED_NONE)
 711       return 0; // not defined !
 712     //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_NONE"));
 713     if (Type==MED_ALL_ELEMENTS)
 714       return _count[_numberOfTypes]-1;
 715     for (med_int i=0; i<_numberOfTypes; i++)
 716       if (_geometricTypes[i]==Type)
 717         return (_count[i+1] - _count[i]);
 718     //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
 719   } else
 720     if (_constituent != NULL)
 721       return _constituent->getNumberOf(Entity,Type);
 722
 723   return 0; // valid if they are nothing else !
 724   //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 725 }
 726
 727 /*! A DOCUMENTER */
 728 //--------------------------------------------------------------//
 729 const med_int* CONNECTIVITY::getValue(medConnectivity TypeConnectivity,
 730                                 medGeometryElement Type)
 731 //--------------------------------------------------------------//
 732 {
 733   if (TypeConnectivity == MED_NODAL)
 734     {
 735           calculateNodalConnectivity();
 736           if (Type==MED_ALL_ELEMENTS)
 737             return _nodal->getValue();
 738           for (med_int i=0; i<_numberOfTypes; i++)
 739             if (_geometricTypes[i]==Type)
 740               return _nodal->getI(_count[i]);
 741     }
 742   else
 743     {
 744       calculateDescendingConnectivity();
 745       if (Type==MED_ALL_ELEMENTS)
 746         return _descending->getValue();
 747       for (med_int i=0; i<_numberOfTypes; i++)
 748         if (_geometricTypes[i]==Type)
 749           return _descending->getI(Type);
 750     }
 751   throw MEDEXCEPTION("Not found");
 752 }
 753
 754 /*! A DOCUMENTER */
 755 //---------------------------------------------------------------------//
 756 const med_int* CONNECTIVITY:: getValueIndex(medConnectivity TypeConnectivity)
 757 //---------------------------------------------------------------------//
 758 {
 759   if (TypeConnectivity == MED_NODAL)
 760     {
 761       calculateNodalConnectivity();
 762       return _nodal->getIndex();
 763     }
 764   else
 765     {
 766       calculateDescendingConnectivity();
 767       return _descending->getIndex();
 768     }
 769 }
 770
 771 /*! Not yet implemented */
 772 //----------------------------------------------//
 773 const med_int* CONNECTIVITY:: getNeighbourhood() const
 774 //----------------------------------------------//
 775 {
 776   throw MEDEXCEPTION("CONNECTIVITY::getNeighbourhood : Not implemented");
 777 }
 778
 779 /*! Returns an array which contains, for each node, all cells
 780     arround it. */
 781 //-------------------------------------------------//
 782 const med_int* CONNECTIVITY::getReverseNodalConnectivity()
 783 //-------------------------------------------------//
 784 {
 785   calculateReverseNodalConnectivity();
 786   return _reverseNodalConnectivity->getValue();
 787 }
 788
 789 /*!  Give index array to use with getReverseConnectivity(MED_NODAL).
 790      It is unusefull with MED_DESCENDING mode, because we have allways two cells.  */
 791 //-------------------------------------------------------//
 792 const med_int* CONNECTIVITY::getReverseNodalConnectivityIndex()
 793 //-------------------------------------------------------//
 794 {
 795   calculateReverseNodalConnectivity();
 796   return _reverseNodalConnectivity->getIndex();
 797 }
 798
 799 /*! Returns an array which contains, for each face (or edge),
 800     the 2 cells of each side. First is cell which face normal is outgoing.
 801     arround it. */
 802 //------------------------------------------------------//
 803 const med_int* CONNECTIVITY::getReverseDescendingConnectivity()
 804 //------------------------------------------------------//
 805 {
 806   // it is in _constituent connectivity only if we are in MED_CELL
 807   // (we could not for instance calculate face-edge connectivity !)
 808   if (_entity!=MED_CELL)
 809     throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivity : Error Only in MED_CELL connectivity");
 810
 811   // we need descending connectivity
 812   calculateDescendingConnectivity();
 813   return _reverseDescendingConnectivity->getValue();
 814 }
 815
 816 /*! calculate the reverse descending Connectivity
 817     and returns the index  ( A DOCUMENTER MIEUX)*/
 818 //-----------------------------------------------------------//
 819 const med_int* CONNECTIVITY::getReverseDescendingConnectivityIndex()
 820 //-----------------------------------------------------------//
 821 {
 822   // it is in _constituent connectivity only if we are in MED_CELL
 823   if (_entity!=MED_CELL)
 824     throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivityIndex : Error Only in MED_CELL connectivity");
 825
 826   // we need descending connectivity
 827   calculateDescendingConnectivity();
 828   return _reverseDescendingConnectivity->getIndex();
 829 }
 830
 831 /*! A DOCUMENTER (et a finir ???) */
 832 //--------------------------------------------//
 833 void CONNECTIVITY::calculateNodalConnectivity()
 834 //--------------------------------------------//
 835 {
 836   if (_nodal==NULL)
 837     {
 838       if (_descending==NULL)
 839         throw MEDEXCEPTION("CONNECTIVITY::calculateNodalConnectivity : No connectivity found !");
 840       // calculate _nodal from _descending
 841     }
 842 }
 843
 844 /*! If not yet done, calculate the nodal Connectivity
 845     and the reverse nodal Connectivity*/
 846 //---------------------------------------------------//
 847 void CONNECTIVITY::calculateReverseNodalConnectivity()
 848 //---------------------------------------------------//
 849 {
 850   const char * LOC = "CONNECTIVITY::calculateReverseNodalConnectivity : ";
 851   BEGIN_OF(LOC);
 852
 853   if (_nodal==NULL)
 854     calculateNodalConnectivity();
 855
 856   if(_reverseNodalConnectivity==NULL) {
 857
 858     med_int node_number = 0;
 859     vector <vector <med_int> > reverse_connectivity;
 860     reverse_connectivity.resize(_numberOfNodes+1);
 861
 862   // Treat all cells types
 863
 864     for (med_int j = 0; j < _numberOfTypes; j++)
 865       {
 866         // node number of the cell type
 867         node_number = _type[j].getNumberOfNodes();
 868         // treat all cells of a particular type
 869         for (med_int k = _count[j]; k < _count[j+1]; k++)
 870           // treat all nodes of the cell type
 871           for (med_int local_node_number = 1; local_node_number < node_number+1; local_node_number++)
 872             {
 873               med_int global_node = _nodal->getIJ(k,local_node_number);
 874               reverse_connectivity[global_node].push_back(k);
 875             }
 876       }
 877
 878     // Full reverse_nodal_connectivity and reverse_nodal_connectivity_index
 879
 880   //calculate size of reverse_nodal_connectivity array
 881     med_int size_reverse_nodal_connectivity  = 0;
 882     for (med_int i = 1; i < _numberOfNodes+1; i++)
 883       size_reverse_nodal_connectivity += reverse_connectivity[i].size();
 884
 885     //MEDSKYLINEARRAY * ReverseConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity);
 886     med_int * reverse_nodal_connectivity_index = new (med_int)[_numberOfNodes+1];
 887     med_int * reverse_nodal_connectivity = new (med_int)[size_reverse_nodal_connectivity];
 888     //const med_int * reverse_nodal_connectivity = ReverseConnectivity->getValue();
 889     //const med_int * reverse_nodal_connectivity_index = ReverseConnectivity->getIndex();
 890
 891     reverse_nodal_connectivity_index[0] = 1;
 892     for (med_int i = 1; i < _numberOfNodes+1; i++)
 893       {
 894         med_int size = reverse_connectivity[i].size();
 895         reverse_nodal_connectivity_index[i] = reverse_nodal_connectivity_index[i-1] + size;
 896         for (med_int j = 0; j < size; j++)
 897           reverse_nodal_connectivity[reverse_nodal_connectivity_index[i-1]-1+j]= reverse_connectivity[i][j];
 898       }
 899
 900     //_reverseNodalConnectivity = ReverseConnectivity;
 901     _reverseNodalConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity,
 902                                                     reverse_nodal_connectivity_index,
 903                                                     reverse_nodal_connectivity);
 904     delete[] reverse_nodal_connectivity_index;
 905     delete[] reverse_nodal_connectivity;
 906  }
 907 }
 908
 909 /*! If not yet done, calculate the Descending Connectivity */
 910 //-------------------------------------------------//
 911 void CONNECTIVITY::calculateDescendingConnectivity()
 912 //-------------------------------------------------//
 913 {
 914   const char * LOC = "CONNECTIVITY::calculateDescendingConnectivity() : ";
 915   BEGIN_OF(LOC);
 916
 917   if (_descending==NULL) {
 918     if (_nodal==NULL) {
 919       MESSAGE(LOC<<"No connectivity found !");
 920       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
 921     }
 922     // calcul _descending from _nodal
 923     // we need CONNECTIVITY for constituent
 924
 925     if (_constituent != NULL)
 926     //      throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"ERROR : No descending connectivity defined, but constituent exist !"));
 927       MESSAGE(LOC<<": No descending connectivity defined, but constituent exist !");
 928
 929     if (_entityDimension == 3)
 930       _constituent = new CONNECTIVITY(MED_FACE);
 931     else if (_entityDimension == 2)
 932       _constituent = new CONNECTIVITY(MED_EDGE);
 933     else {
 934       MESSAGE(LOC<<"We are in 1D");
 935       return;
 936     }
 937     _constituent->_typeConnectivity = MED_NODAL;
 938     _constituent->_numberOfNodes = _numberOfNodes;
 939     // foreach cells, we built array of constituent
 940     int DescendingSize = 0;
 941     for(int i=0; i<_numberOfTypes; i++)
 942       DescendingSize+=(_count[i+1]-_count[i])*_type[i].getNumberOfConstituents(1);
 943     //_descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,DescendingSize);
 944     //const int * descend_connectivity = _descending->getValue();
 945     int * descend_connectivity = new int[DescendingSize];
 946     for (int i=0; i<DescendingSize; i++)
 947       descend_connectivity[i]=0;
 948     //const int * descend_connectivity_index = _descending->getIndex();
 949     int * descend_connectivity_index = new int[_count[_numberOfTypes]];
 950     descend_connectivity_index[0]=1;
 951     medGeometryElement* ConstituentsTypes = new medGeometryElement[2];
 952     ConstituentsTypes[0]=MED_NONE;
 953     ConstituentsTypes[1]=MED_NONE;
 954     int * NumberOfConstituentsForeachType = new int[2];
 955     NumberOfConstituentsForeachType[0]=0;
 956     NumberOfConstituentsForeachType[1]=0;
 957     for(int i=0; i<_numberOfTypes; i++) {
 958       // initialize descend_connectivity_index array :
 959       int NumberOfConstituents = _type[i].getNumberOfConstituents(1);
 960       for (int j=_count[i];j<_count[i+1];j++){
 961         descend_connectivity_index[j]=descend_connectivity_index[j-1]+NumberOfConstituents;
 962         // compute number of constituent of all cell for each type
 963         for(int k=1;k<NumberOfConstituents+1;k++) {
 964           medGeometryElement MEDType = _type[i].getConstituentType(1,k);
 965           if (ConstituentsTypes[0]==MED_NONE) {
 966             ConstituentsTypes[0]=MEDType;
 967             NumberOfConstituentsForeachType[0]++;
 968           } else if (ConstituentsTypes[0]==MEDType)
 969             NumberOfConstituentsForeachType[0]++;
 970           else if (ConstituentsTypes[1]==MED_NONE) {
 971             ConstituentsTypes[1]=MEDType;
 972             NumberOfConstituentsForeachType[1]++;
 973           } else if (ConstituentsTypes[1]==MEDType)
 974             NumberOfConstituentsForeachType[1]++;
 975           else
 976             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<MEDType<<" is different of "<<ConstituentsTypes[0]<<" and "<<ConstituentsTypes[1]));
 977         }
 978       }
 979     }
 980
 981     // we could built _constituent !
 982     int TotalNumberOfConstituents = NumberOfConstituentsForeachType[0]+NumberOfConstituentsForeachType[1];
 983     int TotalNumberOfNodes = NumberOfConstituentsForeachType[0]*(ConstituentsTypes[0]%100)+NumberOfConstituentsForeachType[1]*(ConstituentsTypes[1]%100);
 984
 985     //_constituent->_nodal = new MEDSKYLINEARRAY(TotalNumberOfConstituents,TotalNumberOfNodes);
 986
 987     // we use _constituent->_nodal
 988     //const int * ConstituentNodalConnectivity = _constituent->_nodal->getValue();
 989     //const int * ConstituentNodalConnectivityIndex = _constituent->_nodal->getIndex();
 990     int * ConstituentNodalConnectivity = new int[TotalNumberOfNodes];
 991     int * ConstituentNodalConnectivityIndex = new int[TotalNumberOfConstituents+1];
 992     ConstituentNodalConnectivityIndex[0]=1;
 993
 994     _constituent->_entityDimension=ConstituentsTypes[0]/100;
 995     if (ConstituentsTypes[1]==MED_NONE)
 996       _constituent->_numberOfTypes = 1;
 997     else
 998       _constituent->_numberOfTypes = 2;
 999     _constituent->_geometricTypes = new medGeometryElement[_constituent->_numberOfTypes];
1000     _constituent->_type = new CELLMODEL[_constituent->_numberOfTypes];
1001     _constituent->_count = new med_int[_constituent->_numberOfTypes+1];
1002     _constituent->_count[0]=1;
1003     int* tmp_NumberOfConstituentsForeachType = new med_int[_constituent->_numberOfTypes+1];
1004     tmp_NumberOfConstituentsForeachType[0]=0; // to count constituent of each type
1005     for (int i=0; i<_constituent->_numberOfTypes;i++) {
1006       _constituent->_geometricTypes[i]=ConstituentsTypes[i];
1007       _constituent->_count[i+1]=_constituent->_count[i]+NumberOfConstituentsForeachType[i];
1008       CELLMODEL Type(ConstituentsTypes[i]);
1009       _constituent->_type[i]=Type;
1010       tmp_NumberOfConstituentsForeachType[i+1]=NumberOfConstituentsForeachType[i];
1011       for (int j=tmp_NumberOfConstituentsForeachType[i]; j<tmp_NumberOfConstituentsForeachType[i+1]+tmp_NumberOfConstituentsForeachType[i]; j++)
1012         ConstituentNodalConnectivityIndex[j+1]=ConstituentNodalConnectivityIndex[j]+(ConstituentsTypes[i]%100);
1013     }
1014     delete[] ConstituentsTypes;
1015     delete[] NumberOfConstituentsForeachType;
1016
1017     // we need reverse nodal connectivity
1018     if (! _reverseNodalConnectivity)
1019       calculateReverseNodalConnectivity();
1020     const int * ReverseNodalConnectivityIndex = _reverseNodalConnectivity->getIndex();
1021     const int * ReverseNodalConnectivityValue = _reverseNodalConnectivity->getValue();
1022
1023     // array to keep reverse descending connectivity
1024     int * ReverseDescendingConnectivityValue = new int[TotalNumberOfConstituents*2];
1025
1026     int TotalNumberOfSubCell = 0;
1027     for (int i=0; i<_numberOfTypes; i++) { // loop on all cell type
1028       CELLMODEL Type = _type[i];
1029       int NumberOfNodesPerCell = Type.getNumberOfNodes();
1030       int NumberOfConstituentPerCell = Type.getNumberOfConstituents(1);
1031       for (int j=_count[i]; j<_count[i+1]; j++) // we loop on all cell of this type
1032         for (int k=1; k<=NumberOfConstituentPerCell; k++) { // we loop on all sub cell of it
1033           if (descend_connectivity[descend_connectivity_index[j-1]+k-2]==0) {
1034             // it is a new sub cell !
1035             //      TotalNumberOfSubCell++;
1036             // Which type ?
1037             if (Type.getConstituentType(1,k)==_constituent->_geometricTypes[0]){
1038               tmp_NumberOfConstituentsForeachType[0]++;
1039               TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[0];
1040             } else {
1041               tmp_NumberOfConstituentsForeachType[1]++;
1042               TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[1];
1043             }
1044             //we have maximum two types
1045
1046             descend_connectivity[descend_connectivity_index[j-1]+k-2]=TotalNumberOfSubCell;
1047             ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2]=j;
1048             int NumberOfNodesPerConstituent = Type.getConstituentType(1,k)%100;
1049
1050             int * NodesLists = new int[NumberOfNodesPerConstituent];
1051             for (int l=0; l<NumberOfNodesPerConstituent; l++) {
1052               NodesLists[l]=_nodal->getIJ(j,Type.getNodeConstituent(1,k,l+1));
1053               ConstituentNodalConnectivity[ConstituentNodalConnectivityIndex[TotalNumberOfSubCell-1]-1+l]=NodesLists[l];
1054             }
1055             // we use reverse_nodal_connectivity to find the other element which contain this sub cell
1056
1057             // all elements which contains first node of sub cell :
1058             int ReverseNodalConnectivityIndex_0 = ReverseNodalConnectivityIndex[NodesLists[0]-1];
1059             int ReverseNodalConnectivityIndex_1 = ReverseNodalConnectivityIndex[NodesLists[0]];
1060             int NumberOfCellsInList = ReverseNodalConnectivityIndex_1-ReverseNodalConnectivityIndex_0;
1061
1062             if (NumberOfCellsInList > 0)  { // we could have no element !
1063               int * CellsList = new int[NumberOfCellsInList];
1064               for (int l=ReverseNodalConnectivityIndex_0; l<ReverseNodalConnectivityIndex_1; l++)
1065                 CellsList[l-ReverseNodalConnectivityIndex_0]=ReverseNodalConnectivityValue[l-1];
1066
1067               // foreach node in sub cell, we search elements which are in common
1068               // at the end, we must have only one !
1069
1070               for (int l=1; l<NumberOfNodesPerConstituent; l++) {
1071                 int NewNumberOfCellsInList = 0;
1072                 int * NewCellsList = new int[NumberOfCellsInList];
1073                 for (int l1=0; l1<NumberOfCellsInList; l1++)
1074                   for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<ReverseNodalConnectivityIndex[NodesLists[l]]; l2++) {
1075                     if (CellsList[l1]<ReverseNodalConnectivityValue[l2-1])
1076                       // increasing order : CellsList[l1] are not in elements list of node l
1077                       break;
1078                     if ((CellsList[l1]==ReverseNodalConnectivityValue[l2-1])&(CellsList[l1]!=j)) {
1079                       // we have found one
1080                       NewCellsList[NewNumberOfCellsInList]=CellsList[l1];
1081                       NewNumberOfCellsInList++;
1082                       break;
1083                     }
1084                   }
1085                 NumberOfCellsInList = NewNumberOfCellsInList;
1086
1087                 delete [] CellsList;
1088                 CellsList = NewCellsList;
1089               }
1090
1091               if (NumberOfCellsInList > 0) { // We have found some elements !
1092                 int CellNumber = CellsList[0];
1093                 //delete [] CellsList;
1094                 if (NumberOfCellsInList>1)
1095                   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one other Cell ("<<NumberOfCellsInList<<") !"));
1096
1097                 if (NumberOfCellsInList==1) {
1098                   ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=CellNumber;
1099
1100                   // we search sub cell number in this cell to not calculate it another time
1101                   // which type ?
1102                   CELLMODEL Type2;
1103                   for (int l=0; l<_numberOfTypes; l++)
1104                     if (CellNumber < _count[l+1]) {
1105                       Type2=_type[l];
1106                       break;
1107                     }
1108                   //int sub_cell_count2 = Type2.get_entities_count(1);
1109                   //int nodes_cell_count2 = Type2.get_nodes_count();
1110                   bool find2 = false;
1111                   for (int l=1; l<=Type2.getNumberOfConstituents(1);l++) { // on all sub cell
1112                     int counter = 0;
1113                     for (int m=1; m<=Type2.getConstituentType(1,l)%100; m++)
1114                       for (int n=1; n<=Type.getConstituentType(1,k)%100; n++) {
1115                         if (_nodal->getIJ(CellNumber,Type2.getNodeConstituent(1,l,m)) == NodesLists[n-1])
1116                           counter++;
1117                       }
1118                     if (counter==Type.getConstituentType(1,k)%100) {
1119                       descend_connectivity[descend_connectivity_index[CellNumber-1]+l-2]=-1*TotalNumberOfSubCell; // because, see it in other side !
1120                       find2 = true;
1121                     }
1122                     if (find2)
1123                       break;
1124                   }
1125                   if (!find2)
1126                     INFOS(LOC<<"ERROR ERROR ERROR ERROR ERROR : we find any subcell !!!"); // exception ?
1127                 }
1128               } else {
1129                 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=0;
1130               }
1131               delete[] CellsList;
1132             }
1133
1134             delete[] NodesLists;
1135           }
1136         }
1137     }
1138     // we adjust _constituent
1139     int NumberOfConstituent=0;
1140     int SizeOfConstituentNodal=0;
1141     for (int i=0;i<_constituent->_numberOfTypes; i++) {
1142       NumberOfConstituent += tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1;
1143       SizeOfConstituentNodal += (tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1)*_constituent->_type[i].getNumberOfNodes();
1144     }
1145     // we built new _nodal attribute in _constituent
1146     //MEDSKYLINEARRAY *ConstituentNodal = new MEDSKYLINEARRAY(NumberOfConstituent,SizeOfConstituentNodal);
1147     //const int *ConstituentNodalValue = ConstituentNodal->getValue();
1148     //const int *ConstituentNodalIndex = ConstituentNodal->getIndex();
1149     int *ConstituentNodalValue = new int[SizeOfConstituentNodal];
1150     int *ConstituentNodalIndex = new int[NumberOfConstituent+1];
1151     ConstituentNodalIndex[0]=1;
1152     // we build _reverseDescendingConnectivity
1153     //_reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,2*NumberOfConstituent);
1154     //const int *reverseDescendingConnectivityValue = _reverseDescendingConnectivity->getValue();
1155     //const int *reverseDescendingConnectivityIndex = _reverseDescendingConnectivity->getIndex();
1156     int *reverseDescendingConnectivityValue = new int[2*NumberOfConstituent];
1157     int *reverseDescendingConnectivityIndex = new int[NumberOfConstituent+1];
1158     reverseDescendingConnectivityIndex[0]=1;
1159
1160     // first constituent type
1161     for(int j=0; j<tmp_NumberOfConstituentsForeachType[0]; j++) {
1162       ConstituentNodalIndex[j+1]=ConstituentNodalConnectivityIndex[j+1];
1163       for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1164         ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[k-1];
1165       }
1166       reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1167       for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1168         reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[k-1];
1169       }
1170     }
1171     // second type if any
1172     if (_constituent->_numberOfTypes==2) {
1173       int offset = _constituent->_count[1]-tmp_NumberOfConstituentsForeachType[0]-1;
1174       int offset1=offset*_constituent->_type[0].getNumberOfNodes();
1175       int offset2=offset*2;
1176       int NumberOfNodesPerConstituent = _constituent->_type[1].getNumberOfNodes();
1177       for(int j=tmp_NumberOfConstituentsForeachType[0]; j<NumberOfConstituent; j++) {
1178         ConstituentNodalIndex[j+1]=ConstituentNodalIndex[j]+NumberOfNodesPerConstituent;
1179         for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1180           ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[offset1+k-1];
1181         }
1182         reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1183         for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1184           reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[offset2+k-1];
1185         }
1186       }
1187       _constituent->_count[2]=NumberOfConstituent+1;
1188       // we correct _descending to adjust face number
1189       for(int j=0;j<DescendingSize;j++)
1190         if (descend_connectivity[j]>tmp_NumberOfConstituentsForeachType[0])
1191           descend_connectivity[j]-=offset;
1192
1193     }
1194
1195     delete[] ConstituentNodalConnectivityIndex;
1196     delete[] ConstituentNodalConnectivity;
1197
1198     _descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,
1199                                       DescendingSize,
1200                                       descend_connectivity_index,
1201                                       descend_connectivity);
1202     delete[] descend_connectivity_index;
1203     delete[] descend_connectivity;
1204     _reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,
1205                                                          2*NumberOfConstituent,
1206                                                          reverseDescendingConnectivityIndex,
1207                                                          reverseDescendingConnectivityValue);
1208     delete[] reverseDescendingConnectivityIndex;
1209     delete[] reverseDescendingConnectivityValue;
1210
1211     _constituent->_count[1]=tmp_NumberOfConstituentsForeachType[0]+1;
1212     delete[] tmp_NumberOfConstituentsForeachType;
1213
1214     //delete _constituent->_nodal;
1215     _constituent->_nodal = new MEDSKYLINEARRAY(NumberOfConstituent,
1216                                                SizeOfConstituentNodal,
1217                                                ConstituentNodalIndex,
1218                                                ConstituentNodalValue);
1219
1220     delete[] ConstituentNodalIndex;
1221     delete[] ConstituentNodalValue;
1222
1223     delete[] ReverseDescendingConnectivityValue;
1224   }
1225   END_OF(LOC);
1226 }
1227
1228 //-----------------------------------------------------------------------------------------//
1229 //  void CONNECTIVITY::calculateReverseDescendingConnectivity(CONNECTIVITY *myConnectivity)
1230 //-----------------------------------------------------------------------------------------//
1231 //  {
1232 //    if (_entity==MED_CELL)
1233 //      throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
1234 //        Connectivity : we could not do that with MED_CELL entity !");
1235 //    if (myConnectivity->getEntity()!=MED_CELL)
1236 //      throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
1237 //              Connectivity : we need MED_CELL connectivity !");
1238 //    return;
1239 //  }
1240
1241 /*! Not implemented yet */
1242 //--------------------------------------------------------------------//
1243 void CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity)
1244 //--------------------------------------------------------------------//
1245 {
1246   // Mesh dimension !
1247   return;
1248 }
1249
1250
1251 /*!
1252   Give, for one element number of a specified entity the geometric type
1253   Of this element.
1254
1255   Example : medGeometryElement myType = myConnectivity.getElementType(MED_CELL,35);
1256 */
1257 //--------------------------------------------------------------------//
1258 medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity,int Number) const
1259 //--------------------------------------------------------------------//
1260 {
1261   if (_entity==Entity) {
1262     for(int i=1; i<=_numberOfTypes;i++)
1263       if (Number<_count[i])
1264         return _geometricTypes[i-1];
1265   }
1266   else if (_constituent!=NULL)
1267     return _constituent->getElementType(Entity,Number);
1268   else
1269     throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
1270   throw MEDEXCEPTION(LOCALIZED("getElementType : Wrong Number !"));
1271 }