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     for(int i=0; i<numberOfFamilies; i++) {
 340       FAMILY * myFamily = myFamilies[i] ;
 341       MESSAGE(LOC<<"updating the family (BEGIN) : " << *myFamily);
 342     }
 343
 344     // well we could go !
 345     CONNECTIVITY * oldConstituent = _constituent;
 346
 347 //      for(int i=0; i<numberOfFamilies; i++) {
 348 //        FAMILY * myFamily = myFamilies[i];
 349 //        MESSAGE(LOC<<"updating the family (BEGIN) : " << *myFamily);
 350 //      }
 351
 352     _constituent = (CONNECTIVITY *)NULL;
 353     // for instance we must have nodal connectivity in constituent :
 354     if (oldConstituent->_nodal == NULL)
 355       {
 356         MESSAGE(LOC<<"We have no nodal connectivity of sub cell");
 357         throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
 358       }
 359     int oldNumberOfFace = oldConstituent->_nodal->getNumberOf();
 360     const int * oldConstituentValue = oldConstituent->_nodal->getValue();
 361     const int * oldConstituentIndex = oldConstituent->_nodal->getIndex();
 362
 363     calculateDescendingConnectivity();
 364
 365     int newNumberOfFace = _constituent->_nodal->getNumberOf();
 366     const int * newConstituentValue = _constituent->_nodal->getValue();
 367     const int * newConstituentIndex = _constituent->_nodal->getIndex();
 368
 369     const int * newReverseDescendingIndex =
 370       _reverseDescendingConnectivity->getIndex();
 371
 372     const int * newDescendingIndex = _descending->getIndex();
 373     const int * newDescendingValue = _descending->getValue();
 374
 375     // loop on all family,
 376     //   for all constituent in family, we get it's old connectivity
 377     //   (with oldCconstituentValue and oldConstituentIndex)
 378     //   and search the constituent in newConstituentValue with class
 379     //   ModulusArry
 380     //
 381     //   when a new face is found, replace old constituent
 382     //   number in family with new one
 383     //   If face is not rigth oriented, we must change _descending attribute
 384     //   and _reverseDescendingConnectivity (see calculateDescendingConnectivity()).
 385
 386     // Voila a toi de jouer Nadir :-)
 387
 388     // First we get the renumbering from the oldCconstituentValue and
 389     // oldConstituentIndex in the the new one, newConstituentValue and
 390     // newConstituentIndex with a negative sign if the face is not
 391     // right orented
 392
 393     int * renumberingFromOldToNew = new int [oldNumberOfFace];
 394     int index1 = 0;
 395     int indexm1 = 0;
 396
 397     _constituent->calculateReverseNodalConnectivity();
 398
 399     for (int iOldFace=0;iOldFace<oldNumberOfFace;iOldFace++)
 400       {
 401         int index = 0;
 402
 403         renumberingFromOldToNew[iOldFace] = iOldFace+1;
 404         //        renumberingFromOldToNew[iOldFace] = 999999;
 405
 406         int face_it_beginOld = oldConstituentIndex[iOldFace];
 407         int face_it_endOld = oldConstituentIndex[iOldFace+1];
 408         int face_size_itOld = face_it_endOld - face_it_beginOld;
 409
 410         const int* NodesLists = oldConstituentValue + (face_it_beginOld-1);
 411         int face_size_itNew;
 412
 413         const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity();
 414         const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex();
 415
 416         // set an array wich contains faces numbers arround first node
 417         int BeginIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]-1];
 418         int EndIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]];
 419         int NumberOfFacesInList=EndIndexFaceArrayFirstNode-BeginIndexFaceArrayFirstNode;
 420
 421         int * FacesList = new int[NumberOfFacesInList];
 422
 423         for (int l=BeginIndexFaceArrayFirstNode; l<EndIndexFaceArrayFirstNode; l++){
 424           FacesList[l-BeginIndexFaceArrayFirstNode]=reverseFaceNodal[l-1];
 425         }
 426         // foreach node in sub cell, we search elements which are in common
 427         // at the end, we must have only one !
 428
 429         for (int nodeFaceOld=1; nodeFaceOld<face_size_itOld; nodeFaceOld++)
 430           {
 431             int NewNumberOfFacesInList = 0;
 432             int * NewFacesList = new int[NumberOfFacesInList];
 433
 434             for (int l1=0; l1<NumberOfFacesInList; l1++) {
 435               for (int l2=reverseFaceNodalIndex[NodesLists[nodeFaceOld]-1]; l2<reverseFaceNodalIndex[NodesLists[nodeFaceOld]]; l2++) {
 436                 if (FacesList[l1]<reverseFaceNodal[l2-1])
 437                   // increasing order : FacesList[l1] are not in elements list of node l
 438                   break;
 439                 if (FacesList[l1]==reverseFaceNodal[l2-1]) {
 440                   // we have found one
 441                   NewFacesList[NewNumberOfFacesInList]=FacesList[l1];
 442                   NewNumberOfFacesInList++;
 443                   break;
 444                 }
 445               }
 446             }
 447             NumberOfFacesInList = NewNumberOfFacesInList;
 448             delete [] FacesList;
 449             FacesList = NewFacesList;
 450           }
 451
 452         if (!NumberOfFacesInList==0)
 453           {
 454             if (NumberOfFacesInList>1)
 455               throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one face found ("<<NumberOfFacesInList<<") ! " <<FacesList[0]<<" "<<FacesList[1] ));
 456
 457             MEDMODULUSARRAY modulusArrayOld(face_size_itOld,NodesLists);
 458
 459             int face_it_beginNew = newConstituentIndex[FacesList[0]-1];
 460             int face_it_endNew = newConstituentIndex[FacesList[0]];
 461             face_size_itNew = face_it_endNew - face_it_beginNew;
 462
 463             const int * newNodesLists = newConstituentValue+newConstituentIndex[FacesList[0]-1]-1;
 464             MEDMODULUSARRAY modulusArrayNew(face_size_itNew,newNodesLists);
 465
 466             int retCompareNewOld = modulusArrayNew.compare(modulusArrayOld);
 467
 468             //SCRUTE(retCompareNewOld);
 469
 470             // Real new face found
 471
 472             if(retCompareNewOld == 1)
 473               {
 474                 renumberingFromOldToNew[iOldFace] = FacesList[0];
 475                 index = 1;
 476                 index1++;
 477               }
 478
 479             // Reverse new face found
 480
 481             if(retCompareNewOld == -1)
 482               {
 483                 renumberingFromOldToNew[iOldFace] = FacesList[0];
 484                 index = 1;
 485                 indexm1++;
 486
 487                 int face_it_begin = newReverseDescendingIndex[FacesList[0]-1];
 488                 int face_it_end = newReverseDescendingIndex[FacesList[0]];
 489                 int face_size_it = face_it_end - face_it_begin;
 490
 491                 if (face_size_it == 1)
 492                   throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
 493
 494                 if (face_size_it > 2)
 495                   throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This face/edge should not be a (d-1) cell because it has "<<face_size_it<<" neighbouring"));
 496
 497                 // we have always 2 neighbourings
 498                 int cell1 = _reverseDescendingConnectivity->getIJ(FacesList[0],1);
 499                 int cell2 = _reverseDescendingConnectivity->getIJ(FacesList[0],2);
 500                 // PROVISOIRE : en attendant que le SKYLINEARRAY de ReverseDescending soit correct (sans le zero)
 501                 //                  if (cell2 == 0)
 502                 //                    throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
 503
 504                 if (cell2 != 0) { // we are not on border !!!!
 505
 506                   _reverseDescendingConnectivity->setIJ(FacesList[0],1,cell2);
 507                   // Updating _constituent->_nodal because of reversity
 508                   const int * oldArray = oldConstituentValue+face_it_beginOld-1;
 509                   for(int iarray=1;iarray<=face_size_itNew;iarray++){
 510                     _constituent->_nodal->setIJ(FacesList[0],iarray,oldArray[iarray-1]);
 511                   }
 512
 513                   // Updating _reverseDescendingConnectivity
 514
 515
 516                   _reverseDescendingConnectivity->setIJ(FacesList[0],2,cell1);
 517
 518                   // Updating _descending for cell1 and cell2
 519                   for(int iface=newDescendingIndex[cell1-1];iface<=newDescendingIndex[cell1];iface++)
 520                     if (_descending->getIndexValue(iface)==FacesList[0])
 521                       _descending->setIndexValue(iface,-FacesList[0]);
 522                     else if (_descending->getIndexValue(iface)==-FacesList[0])
 523                       _descending->setIndexValue(iface,FacesList[0]);
 524                   // else nothing to do
 525
 526                   for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];iface++)
 527                     if (_descending->getIndexValue(iface)==FacesList[0])
 528                       _descending->setIndexValue(iface,-FacesList[0]);
 529                     else if (_descending->getIndexValue(iface)==-FacesList[0])
 530                       _descending->setIndexValue(iface,FacesList[0]);
 531                   // else nothing to do
 532
 533                 } else {// else we are on border and we do nothing !!!!!!!!
 534                   INFOS("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
 535                   INFOS(LOC<<" Boudary FACE "<<iOldFace+1<<" are wrong oriented !");
 536                   INFOS("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
 537                 }
 538               }
 539
 540             if(index == 0)
 541               {
 542                 INFOS(LOC<<"Renumbering problem with the Face connectivity given by the User and the new Connectivity computed");
 543                 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a Face connectivity problem"));
 544               }
 545           }
 546         delete[] FacesList;
 547       }
 548
 549     MESSAGE(LOC<<"The Renumbering is finished and the status is");
 550
 551     // Updating the Family
 552
 553     for(int i=0; i<numberOfFamilies; i++) {
 554       FAMILY * myFamily = myFamilies[i];
 555
 556       MEDSKYLINEARRAY * number = myFamily->getnumber();
 557       int numberOfLines_skyline = number->getNumberOf();
 558       const int * index_skyline = number->getIndex();
 559
 560       for (int i=0;i<numberOfLines_skyline;i++) {
 561         for (int j=index_skyline[i]; j<index_skyline[i+1];j++) {
 562           number->setIndexValue(j,renumberingFromOldToNew[number->getIndexValue(j)-1]);
 563         }
 564       }
 565       MESSAGE(LOC<<"updating the family (END) : " << *myFamily);
 566     }
 567
 568     delete oldConstituent ;
 569     delete [] renumberingFromOldToNew;
 570   }
 571
 572
 573   END_OF(LOC);
 574   return;
 575 }
 576
 577 //------------------------------------------------------------------------------------------------------------------//
 578 const med_int * CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type)
 579 //------------------------------------------------------------------------------------------------------------------//
 580 {
 581   const char * LOC = "CONNECTIVITY::getConnectivity";
 582   BEGIN_OF(LOC);
 583
 584   MEDSKYLINEARRAY * Connectivity;
 585   if (Entity==_entity) {
 586
 587     if (ConnectivityType==MED_NODAL)
 588       {
 589         calculateNodalConnectivity();
 590         Connectivity=_nodal;
 591       }
 592     else
 593       {
 594         calculateDescendingConnectivity();
 595         Connectivity=_descending;
 596       }
 597
 598     if (Connectivity!=NULL)
 599       if (Type==MED_ALL_ELEMENTS)
 600         return Connectivity->getValue();
 601       else {
 602         for (med_int i=0; i<_numberOfTypes; i++)
 603           if (_geometricTypes[i]==Type)
 604             //return Connectivity->getI(i+1);
 605             return Connectivity->getI(_count[i]);
 606         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
 607       }
 608     else
 609       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
 610   } else
 611     if (_constituent != NULL)
 612       return _constituent->getConnectivity(ConnectivityType,Entity,Type);
 613
 614   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 615 }
 616
 617 /*!  Give morse index array to use with
 618      getConnectivity(MED_FULL_INTERLACE,mode,entity,MED_ALL_ELEMENTS).
 619
 620       Each value give start index for corresponding entity in connectivity array.
 621
 622       Example : i-th element, j-th node of it :
 623       - In C mode : Connectivity[ConnectivityIndex[i]-1+j-1]
 624       - In fortran mode : Connectivity[ConnectivityIndex[i]+j] */
 625 //-----------------------------------------------------------------------------------------------//
 626 const med_int * CONNECTIVITY::getConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity)
 627 //------------------------------------------------------------------------------------------------//
 628 {
 629   const char * LOC = "CONNECTIVITY::getConnectivityIndex";
 630   BEGIN_OF(LOC);
 631
 632   MEDSKYLINEARRAY * Connectivity;
 633   if (Entity==_entity) {
 634
 635     if (ConnectivityType==MED_NODAL)
 636       Connectivity=_nodal;
 637     else
 638       Connectivity=_descending;
 639
 640     if (Connectivity!=NULL)
 641       return Connectivity->getIndex();
 642     else
 643       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
 644
 645   } else
 646     if (_constituent != NULL)
 647       return _constituent->getConnectivityIndex(ConnectivityType,Entity);
 648
 649   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 650 }
 651
 652 /*! A DOCUMENTER */
 653 //--------------------------------------------------------------//
 654 const CELLMODEL & CONNECTIVITY::getType(medGeometryElement Type) const
 655 //--------------------------------------------------------------//
 656 {
 657   const char * LOC = "CONNECTIVITY::getType";
 658   BEGIN_OF(LOC);
 659
 660   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
 661     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE !"));
 662   for (med_int i=0; i<_numberOfTypes; i++)
 663     if (_geometricTypes[i]==Type)
 664       return _type[i];
 665   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
 666 }
 667
 668 /*!  Returns the number of elements of type <med geometrie element>.
 669      Note : not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
 670 //------------------------------------------------------------------------//
 671 med_int CONNECTIVITY::getNumberOfNodesInType(medGeometryElement Type) const
 672 //------------------------------------------------------------------------//
 673 {
 674   const char * LOC = "CONNECTIVITY::getNumberOfNodesInType";
 675   BEGIN_OF(LOC);
 676
 677   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
 678     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!"));
 679   for (med_int i=0; i<_numberOfTypes; i++)
 680     if (_geometricTypes[i]==Type)
 681       return _type[i].getNumberOfNodes();
 682   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
 683 }
 684
 685 /*!  Returns the number of geometric sub cells of <med geometrie element> type.
 686 not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
 687 //------------------------------------------------------------------------//
 688 med_int CONNECTIVITY::getNumberOfSubCellInType(medGeometryElement Type) const
 689 //------------------------------------------------------------------------//
 690 {
 691   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
 692     throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!");
 693   for (med_int i=0; i<_numberOfTypes; i++)
 694     if (_geometricTypes[i]==Type)
 695       return _type[i].getNumberOfConstituents(1);
 696   throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement not found !");
 697 }
 698
 699 /*!  Returns the number of elements of type <med geometrie element>.
 700
 701      Note :
 702       - Implemented for MED_ALL_ELEMENTS
 703       - Not implemented for MED_ALL_ENTITIES (A VERIFIER)
 704       - Not implemented for MED_NONE */
 705 //-----------------------------------------------------------------------------------//
 706 med_int CONNECTIVITY::getNumberOf(medEntityMesh Entity, medGeometryElement Type) const
 707 //-----------------------------------------------------------------------------------//
 708 {
 709   const char * LOC = "CONNECTIVITY::getNumberOf";
 710   BEGIN_OF(LOC);
 711
 712   MESSAGE(LOC<<" Entity = "<< Entity << ", _entity = "<<_entity);
 713
 714   if (Entity==_entity) {
 715     if (Type==MED_NONE)
 716       return 0; // not defined !
 717     //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_NONE"));
 718     if (Type==MED_ALL_ELEMENTS)
 719       return _count[_numberOfTypes]-1;
 720     for (med_int i=0; i<_numberOfTypes; i++)
 721       if (_geometricTypes[i]==Type)
 722         return (_count[i+1] - _count[i]);
 723     //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
 724   } else
 725     if (_constituent != NULL)
 726       return _constituent->getNumberOf(Entity,Type);
 727
 728   return 0; // valid if they are nothing else !
 729   //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 730 }
 731
 732 /*! A DOCUMENTER */
 733 //--------------------------------------------------------------//
 734 const med_int* CONNECTIVITY::getValue(medConnectivity TypeConnectivity,
 735                                 medGeometryElement Type)
 736 //--------------------------------------------------------------//
 737 {
 738   if (TypeConnectivity == MED_NODAL)
 739     {
 740           calculateNodalConnectivity();
 741           if (Type==MED_ALL_ELEMENTS)
 742             return _nodal->getValue();
 743           for (med_int i=0; i<_numberOfTypes; i++)
 744             if (_geometricTypes[i]==Type)
 745               return _nodal->getI(_count[i]);
 746     }
 747   else
 748     {
 749       calculateDescendingConnectivity();
 750       if (Type==MED_ALL_ELEMENTS)
 751         return _descending->getValue();
 752       for (med_int i=0; i<_numberOfTypes; i++)
 753         if (_geometricTypes[i]==Type)
 754           return _descending->getI(Type);
 755     }
 756   throw MEDEXCEPTION("Not found");
 757 }
 758
 759 /*! A DOCUMENTER */
 760 //---------------------------------------------------------------------//
 761 const med_int* CONNECTIVITY:: getValueIndex(medConnectivity TypeConnectivity)
 762 //---------------------------------------------------------------------//
 763 {
 764   if (TypeConnectivity == MED_NODAL)
 765     {
 766       calculateNodalConnectivity();
 767       return _nodal->getIndex();
 768     }
 769   else
 770     {
 771       calculateDescendingConnectivity();
 772       return _descending->getIndex();
 773     }
 774 }
 775
 776 /*! Not yet implemented */
 777 //----------------------------------------------//
 778 const med_int* CONNECTIVITY:: getNeighbourhood() const
 779 //----------------------------------------------//
 780 {
 781   throw MEDEXCEPTION("CONNECTIVITY::getNeighbourhood : Not implemented");
 782 }
 783
 784 /*! Returns an array which contains, for each node, all cells
 785     arround it. */
 786 //-------------------------------------------------//
 787 const med_int* CONNECTIVITY::getReverseNodalConnectivity()
 788 //-------------------------------------------------//
 789 {
 790   calculateReverseNodalConnectivity();
 791   return _reverseNodalConnectivity->getValue();
 792 }
 793
 794 /*!  Give index array to use with getReverseConnectivity(MED_NODAL).
 795      It is unusefull with MED_DESCENDING mode, because we have allways two cells.  */
 796 //-------------------------------------------------------//
 797 const med_int* CONNECTIVITY::getReverseNodalConnectivityIndex()
 798 //-------------------------------------------------------//
 799 {
 800   calculateReverseNodalConnectivity();
 801   return _reverseNodalConnectivity->getIndex();
 802 }
 803
 804 /*! Returns an array which contains, for each face (or edge),
 805     the 2 cells of each side. First is cell which face normal is outgoing.
 806     arround it. */
 807 //------------------------------------------------------//
 808 const med_int* CONNECTIVITY::getReverseDescendingConnectivity()
 809 //------------------------------------------------------//
 810 {
 811   // it is in _constituent connectivity only if we are in MED_CELL
 812   // (we could not for instance calculate face-edge connectivity !)
 813   if (_entity!=MED_CELL)
 814     throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivity : Error Only in MED_CELL connectivity");
 815
 816   // we need descending connectivity
 817   calculateDescendingConnectivity();
 818   return _reverseDescendingConnectivity->getValue();
 819 }
 820
 821 /*! calculate the reverse descending Connectivity
 822     and returns the index  ( A DOCUMENTER MIEUX)*/
 823 //-----------------------------------------------------------//
 824 const med_int* CONNECTIVITY::getReverseDescendingConnectivityIndex()
 825 //-----------------------------------------------------------//
 826 {
 827   // it is in _constituent connectivity only if we are in MED_CELL
 828   if (_entity!=MED_CELL)
 829     throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivityIndex : Error Only in MED_CELL connectivity");
 830
 831   // we need descending connectivity
 832   calculateDescendingConnectivity();
 833   return _reverseDescendingConnectivity->getIndex();
 834 }
 835
 836 /*! A DOCUMENTER (et a finir ???) */
 837 //--------------------------------------------//
 838 void CONNECTIVITY::calculateNodalConnectivity()
 839 //--------------------------------------------//
 840 {
 841   if (_nodal==NULL)
 842     {
 843       if (_descending==NULL)
 844         throw MEDEXCEPTION("CONNECTIVITY::calculateNodalConnectivity : No connectivity found !");
 845       // calculate _nodal from _descending
 846     }
 847 }
 848
 849 /*! If not yet done, calculate the nodal Connectivity
 850     and the reverse nodal Connectivity*/
 851 //---------------------------------------------------//
 852 void CONNECTIVITY::calculateReverseNodalConnectivity()
 853 //---------------------------------------------------//
 854 {
 855   const char * LOC = "CONNECTIVITY::calculateReverseNodalConnectivity : ";
 856   BEGIN_OF(LOC);
 857
 858   if (_nodal==NULL)
 859     calculateNodalConnectivity();
 860
 861   if(_reverseNodalConnectivity==NULL) {
 862
 863     med_int node_number = 0;
 864     vector <vector <med_int> > reverse_connectivity;
 865     reverse_connectivity.resize(_numberOfNodes+1);
 866
 867   // Treat all cells types
 868
 869     for (med_int j = 0; j < _numberOfTypes; j++)
 870       {
 871         // node number of the cell type
 872         node_number = _type[j].getNumberOfNodes();
 873         // treat all cells of a particular type
 874         for (med_int k = _count[j]; k < _count[j+1]; k++)
 875           // treat all nodes of the cell type
 876           for (med_int local_node_number = 1; local_node_number < node_number+1; local_node_number++)
 877             {
 878               med_int global_node = _nodal->getIJ(k,local_node_number);
 879               reverse_connectivity[global_node].push_back(k);
 880             }
 881       }
 882
 883     // Full reverse_nodal_connectivity and reverse_nodal_connectivity_index
 884
 885   //calculate size of reverse_nodal_connectivity array
 886     med_int size_reverse_nodal_connectivity  = 0;
 887     for (med_int i = 1; i < _numberOfNodes+1; i++)
 888       size_reverse_nodal_connectivity += reverse_connectivity[i].size();
 889
 890     //MEDSKYLINEARRAY * ReverseConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity);
 891     med_int * reverse_nodal_connectivity_index = new (med_int)[_numberOfNodes+1];
 892     med_int * reverse_nodal_connectivity = new (med_int)[size_reverse_nodal_connectivity];
 893     //const med_int * reverse_nodal_connectivity = ReverseConnectivity->getValue();
 894     //const med_int * reverse_nodal_connectivity_index = ReverseConnectivity->getIndex();
 895
 896     reverse_nodal_connectivity_index[0] = 1;
 897     for (med_int i = 1; i < _numberOfNodes+1; i++)
 898       {
 899         med_int size = reverse_connectivity[i].size();
 900         reverse_nodal_connectivity_index[i] = reverse_nodal_connectivity_index[i-1] + size;
 901         for (med_int j = 0; j < size; j++)
 902           reverse_nodal_connectivity[reverse_nodal_connectivity_index[i-1]-1+j]= reverse_connectivity[i][j];
 903       }
 904
 905     //_reverseNodalConnectivity = ReverseConnectivity;
 906     _reverseNodalConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity,
 907                                                     reverse_nodal_connectivity_index,
 908                                                     reverse_nodal_connectivity);
 909     delete[] reverse_nodal_connectivity_index;
 910     delete[] reverse_nodal_connectivity;
 911  }
 912 }
 913
 914 /*! If not yet done, calculate the Descending Connectivity */
 915 //-------------------------------------------------//
 916 void CONNECTIVITY::calculateDescendingConnectivity()
 917 //-------------------------------------------------//
 918 {
 919   const char * LOC = "CONNECTIVITY::calculateDescendingConnectivity() : ";
 920   BEGIN_OF(LOC);
 921
 922   if (_descending==NULL) {
 923     if (_nodal==NULL) {
 924       MESSAGE(LOC<<"No connectivity found !");
 925       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
 926     }
 927     // calcul _descending from _nodal
 928     // we need CONNECTIVITY for constituent
 929
 930     if (_constituent != NULL)
 931     //      throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"ERROR : No descending connectivity defined, but constituent exist !"));
 932       MESSAGE(LOC<<": No descending connectivity defined, but constituent exist !");
 933
 934     if (_entityDimension == 3)
 935       _constituent = new CONNECTIVITY(MED_FACE);
 936     else if (_entityDimension == 2)
 937       _constituent = new CONNECTIVITY(MED_EDGE);
 938     else {
 939       MESSAGE(LOC<<"We are in 1D");
 940       return;
 941     }
 942     _constituent->_typeConnectivity = MED_NODAL;
 943     _constituent->_numberOfNodes = _numberOfNodes;
 944     // foreach cells, we built array of constituent
 945     int DescendingSize = 0;
 946     for(int i=0; i<_numberOfTypes; i++)
 947       DescendingSize+=(_count[i+1]-_count[i])*_type[i].getNumberOfConstituents(1);
 948     //_descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,DescendingSize);
 949     //const int * descend_connectivity = _descending->getValue();
 950     int * descend_connectivity = new int[DescendingSize];
 951     for (int i=0; i<DescendingSize; i++)
 952       descend_connectivity[i]=0;
 953     //const int * descend_connectivity_index = _descending->getIndex();
 954     int * descend_connectivity_index = new int[_count[_numberOfTypes]];
 955     descend_connectivity_index[0]=1;
 956     medGeometryElement* ConstituentsTypes = new medGeometryElement[2];
 957     ConstituentsTypes[0]=MED_NONE;
 958     ConstituentsTypes[1]=MED_NONE;
 959     int * NumberOfConstituentsForeachType = new int[2];
 960     NumberOfConstituentsForeachType[0]=0;
 961     NumberOfConstituentsForeachType[1]=0;
 962     for(int i=0; i<_numberOfTypes; i++) {
 963       // initialize descend_connectivity_index array :
 964       int NumberOfConstituents = _type[i].getNumberOfConstituents(1);
 965       for (int j=_count[i];j<_count[i+1];j++){
 966         descend_connectivity_index[j]=descend_connectivity_index[j-1]+NumberOfConstituents;
 967         // compute number of constituent of all cell for each type
 968         for(int k=1;k<NumberOfConstituents+1;k++) {
 969           medGeometryElement MEDType = _type[i].getConstituentType(1,k);
 970           if (ConstituentsTypes[0]==MED_NONE) {
 971             ConstituentsTypes[0]=MEDType;
 972             NumberOfConstituentsForeachType[0]++;
 973           } else if (ConstituentsTypes[0]==MEDType)
 974             NumberOfConstituentsForeachType[0]++;
 975           else if (ConstituentsTypes[1]==MED_NONE) {
 976             ConstituentsTypes[1]=MEDType;
 977             NumberOfConstituentsForeachType[1]++;
 978           } else if (ConstituentsTypes[1]==MEDType)
 979             NumberOfConstituentsForeachType[1]++;
 980           else
 981             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<MEDType<<" is different of "<<ConstituentsTypes[0]<<" and "<<ConstituentsTypes[1]));
 982         }
 983       }
 984     }
 985
 986     // we could built _constituent !
 987     int TotalNumberOfConstituents = NumberOfConstituentsForeachType[0]+NumberOfConstituentsForeachType[1];
 988     int TotalNumberOfNodes = NumberOfConstituentsForeachType[0]*(ConstituentsTypes[0]%100)+NumberOfConstituentsForeachType[1]*(ConstituentsTypes[1]%100);
 989
 990     //_constituent->_nodal = new MEDSKYLINEARRAY(TotalNumberOfConstituents,TotalNumberOfNodes);
 991
 992     // we use _constituent->_nodal
 993     //const int * ConstituentNodalConnectivity = _constituent->_nodal->getValue();
 994     //const int * ConstituentNodalConnectivityIndex = _constituent->_nodal->getIndex();
 995     int * ConstituentNodalConnectivity = new int[TotalNumberOfNodes];
 996     int * ConstituentNodalConnectivityIndex = new int[TotalNumberOfConstituents+1];
 997     ConstituentNodalConnectivityIndex[0]=1;
 998
 999     _constituent->_entityDimension=ConstituentsTypes[0]/100;
1000     if (ConstituentsTypes[1]==MED_NONE)
1001       _constituent->_numberOfTypes = 1;
1002     else
1003       _constituent->_numberOfTypes = 2;
1004     _constituent->_geometricTypes = new medGeometryElement[_constituent->_numberOfTypes];
1005     _constituent->_type = new CELLMODEL[_constituent->_numberOfTypes];
1006     _constituent->_count = new med_int[_constituent->_numberOfTypes+1];
1007     _constituent->_count[0]=1;
1008     int* tmp_NumberOfConstituentsForeachType = new med_int[_constituent->_numberOfTypes+1];
1009     tmp_NumberOfConstituentsForeachType[0]=0; // to count constituent of each type
1010     for (int i=0; i<_constituent->_numberOfTypes;i++) {
1011       _constituent->_geometricTypes[i]=ConstituentsTypes[i];
1012       _constituent->_count[i+1]=_constituent->_count[i]+NumberOfConstituentsForeachType[i];
1013       CELLMODEL Type(ConstituentsTypes[i]);
1014       _constituent->_type[i]=Type;
1015       tmp_NumberOfConstituentsForeachType[i+1]=NumberOfConstituentsForeachType[i];
1016       for (int j=tmp_NumberOfConstituentsForeachType[i]; j<tmp_NumberOfConstituentsForeachType[i+1]+tmp_NumberOfConstituentsForeachType[i]; j++)
1017         ConstituentNodalConnectivityIndex[j+1]=ConstituentNodalConnectivityIndex[j]+(ConstituentsTypes[i]%100);
1018     }
1019     delete[] ConstituentsTypes;
1020     delete[] NumberOfConstituentsForeachType;
1021
1022     // we need reverse nodal connectivity
1023     if (! _reverseNodalConnectivity)
1024       calculateReverseNodalConnectivity();
1025     const int * ReverseNodalConnectivityIndex = _reverseNodalConnectivity->getIndex();
1026     const int * ReverseNodalConnectivityValue = _reverseNodalConnectivity->getValue();
1027
1028     // array to keep reverse descending connectivity
1029     int * ReverseDescendingConnectivityValue = new int[TotalNumberOfConstituents*2];
1030
1031     int TotalNumberOfSubCell = 0;
1032     for (int i=0; i<_numberOfTypes; i++) { // loop on all cell type
1033       CELLMODEL Type = _type[i];
1034       int NumberOfNodesPerCell = Type.getNumberOfNodes();
1035       int NumberOfConstituentPerCell = Type.getNumberOfConstituents(1);
1036       for (int j=_count[i]; j<_count[i+1]; j++) // we loop on all cell of this type
1037         for (int k=1; k<=NumberOfConstituentPerCell; k++) { // we loop on all sub cell of it
1038           if (descend_connectivity[descend_connectivity_index[j-1]+k-2]==0) {
1039             // it is a new sub cell !
1040             //      TotalNumberOfSubCell++;
1041             // Which type ?
1042             if (Type.getConstituentType(1,k)==_constituent->_geometricTypes[0]){
1043               tmp_NumberOfConstituentsForeachType[0]++;
1044               TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[0];
1045             } else {
1046               tmp_NumberOfConstituentsForeachType[1]++;
1047               TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[1];
1048             }
1049             //we have maximum two types
1050
1051             descend_connectivity[descend_connectivity_index[j-1]+k-2]=TotalNumberOfSubCell;
1052             ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2]=j;
1053             int NumberOfNodesPerConstituent = Type.getConstituentType(1,k)%100;
1054
1055             int * NodesLists = new int[NumberOfNodesPerConstituent];
1056             for (int l=0; l<NumberOfNodesPerConstituent; l++) {
1057               NodesLists[l]=_nodal->getIJ(j,Type.getNodeConstituent(1,k,l+1));
1058               ConstituentNodalConnectivity[ConstituentNodalConnectivityIndex[TotalNumberOfSubCell-1]-1+l]=NodesLists[l];
1059             }
1060             // we use reverse_nodal_connectivity to find the other element which contain this sub cell
1061
1062             // all elements which contains first node of sub cell :
1063             int ReverseNodalConnectivityIndex_0 = ReverseNodalConnectivityIndex[NodesLists[0]-1];
1064             int ReverseNodalConnectivityIndex_1 = ReverseNodalConnectivityIndex[NodesLists[0]];
1065             int NumberOfCellsInList = ReverseNodalConnectivityIndex_1-ReverseNodalConnectivityIndex_0;
1066
1067             if (NumberOfCellsInList > 0)  { // we could have no element !
1068               int * CellsList = new int[NumberOfCellsInList];
1069               for (int l=ReverseNodalConnectivityIndex_0; l<ReverseNodalConnectivityIndex_1; l++)
1070                 CellsList[l-ReverseNodalConnectivityIndex_0]=ReverseNodalConnectivityValue[l-1];
1071
1072               // foreach node in sub cell, we search elements which are in common
1073               // at the end, we must have only one !
1074
1075               for (int l=1; l<NumberOfNodesPerConstituent; l++) {
1076                 int NewNumberOfCellsInList = 0;
1077                 int * NewCellsList = new int[NumberOfCellsInList];
1078                 for (int l1=0; l1<NumberOfCellsInList; l1++)
1079                   for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<ReverseNodalConnectivityIndex[NodesLists[l]]; l2++) {
1080                     if (CellsList[l1]<ReverseNodalConnectivityValue[l2-1])
1081                       // increasing order : CellsList[l1] are not in elements list of node l
1082                       break;
1083                     if ((CellsList[l1]==ReverseNodalConnectivityValue[l2-1])&(CellsList[l1]!=j)) {
1084                       // we have found one
1085                       NewCellsList[NewNumberOfCellsInList]=CellsList[l1];
1086                       NewNumberOfCellsInList++;
1087                       break;
1088                     }
1089                   }
1090                 NumberOfCellsInList = NewNumberOfCellsInList;
1091
1092                 delete [] CellsList;
1093                 CellsList = NewCellsList;
1094               }
1095
1096               if (NumberOfCellsInList > 0) { // We have found some elements !
1097                 int CellNumber = CellsList[0];
1098                 //delete [] CellsList;
1099                 if (NumberOfCellsInList>1)
1100                   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one other Cell ("<<NumberOfCellsInList<<") !"));
1101
1102                 if (NumberOfCellsInList==1) {
1103                   ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=CellNumber;
1104
1105                   // we search sub cell number in this cell to not calculate it another time
1106                   // which type ?
1107                   CELLMODEL Type2;
1108                   for (int l=0; l<_numberOfTypes; l++)
1109                     if (CellNumber < _count[l+1]) {
1110                       Type2=_type[l];
1111                       break;
1112                     }
1113                   //int sub_cell_count2 = Type2.get_entities_count(1);
1114                   //int nodes_cell_count2 = Type2.get_nodes_count();
1115                   bool find2 = false;
1116                   for (int l=1; l<=Type2.getNumberOfConstituents(1);l++) { // on all sub cell
1117                     int counter = 0;
1118                     for (int m=1; m<=Type2.getConstituentType(1,l)%100; m++)
1119                       for (int n=1; n<=Type.getConstituentType(1,k)%100; n++) {
1120                         if (_nodal->getIJ(CellNumber,Type2.getNodeConstituent(1,l,m)) == NodesLists[n-1])
1121                           counter++;
1122                       }
1123                     if (counter==Type.getConstituentType(1,k)%100) {
1124                       descend_connectivity[descend_connectivity_index[CellNumber-1]+l-2]=-1*TotalNumberOfSubCell; // because, see it in other side !
1125                       find2 = true;
1126                     }
1127                     if (find2)
1128                       break;
1129                   }
1130                   if (!find2)
1131                     INFOS(LOC<<"ERROR ERROR ERROR ERROR ERROR : we find any subcell !!!"); // exception ?
1132                 }
1133               } else {
1134                 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=0;
1135               }
1136               delete[] CellsList;
1137             }
1138
1139             delete[] NodesLists;
1140           }
1141         }
1142     }
1143     // we adjust _constituent
1144     int NumberOfConstituent=0;
1145     int SizeOfConstituentNodal=0;
1146     for (int i=0;i<_constituent->_numberOfTypes; i++) {
1147       NumberOfConstituent += tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1;
1148       SizeOfConstituentNodal += (tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1)*_constituent->_type[i].getNumberOfNodes();
1149     }
1150     // we built new _nodal attribute in _constituent
1151     //MEDSKYLINEARRAY *ConstituentNodal = new MEDSKYLINEARRAY(NumberOfConstituent,SizeOfConstituentNodal);
1152     //const int *ConstituentNodalValue = ConstituentNodal->getValue();
1153     //const int *ConstituentNodalIndex = ConstituentNodal->getIndex();
1154     int *ConstituentNodalValue = new int[SizeOfConstituentNodal];
1155     int *ConstituentNodalIndex = new int[NumberOfConstituent+1];
1156     ConstituentNodalIndex[0]=1;
1157     // we build _reverseDescendingConnectivity
1158     //_reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,2*NumberOfConstituent);
1159     //const int *reverseDescendingConnectivityValue = _reverseDescendingConnectivity->getValue();
1160     //const int *reverseDescendingConnectivityIndex = _reverseDescendingConnectivity->getIndex();
1161     int *reverseDescendingConnectivityValue = new int[2*NumberOfConstituent];
1162     int *reverseDescendingConnectivityIndex = new int[NumberOfConstituent+1];
1163     reverseDescendingConnectivityIndex[0]=1;
1164
1165     // first constituent type
1166     for(int j=0; j<tmp_NumberOfConstituentsForeachType[0]; j++) {
1167       ConstituentNodalIndex[j+1]=ConstituentNodalConnectivityIndex[j+1];
1168       for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1169         ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[k-1];
1170       }
1171       reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1172       for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1173         reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[k-1];
1174       }
1175     }
1176     // second type if any
1177     if (_constituent->_numberOfTypes==2) {
1178       int offset = _constituent->_count[1]-tmp_NumberOfConstituentsForeachType[0]-1;
1179       int offset1=offset*_constituent->_type[0].getNumberOfNodes();
1180       int offset2=offset*2;
1181       int NumberOfNodesPerConstituent = _constituent->_type[1].getNumberOfNodes();
1182       for(int j=tmp_NumberOfConstituentsForeachType[0]; j<NumberOfConstituent; j++) {
1183         ConstituentNodalIndex[j+1]=ConstituentNodalIndex[j]+NumberOfNodesPerConstituent;
1184         for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1185           ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[offset1+k-1];
1186         }
1187         reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1188         for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1189           reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[offset2+k-1];
1190         }
1191       }
1192       _constituent->_count[2]=NumberOfConstituent+1;
1193       // we correct _descending to adjust face number
1194       for(int j=0;j<DescendingSize;j++)
1195         if (descend_connectivity[j]>tmp_NumberOfConstituentsForeachType[0])
1196           descend_connectivity[j]-=offset;
1197
1198     }
1199
1200     delete[] ConstituentNodalConnectivityIndex;
1201     delete[] ConstituentNodalConnectivity;
1202
1203     _descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,
1204                                       DescendingSize,
1205                                       descend_connectivity_index,
1206                                       descend_connectivity);
1207     delete[] descend_connectivity_index;
1208     delete[] descend_connectivity;
1209     _reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,
1210                                                          2*NumberOfConstituent,
1211                                                          reverseDescendingConnectivityIndex,
1212                                                          reverseDescendingConnectivityValue);
1213     delete[] reverseDescendingConnectivityIndex;
1214     delete[] reverseDescendingConnectivityValue;
1215
1216     _constituent->_count[1]=tmp_NumberOfConstituentsForeachType[0]+1;
1217     delete[] tmp_NumberOfConstituentsForeachType;
1218
1219     //delete _constituent->_nodal;
1220     _constituent->_nodal = new MEDSKYLINEARRAY(NumberOfConstituent,
1221                                                SizeOfConstituentNodal,
1222                                                ConstituentNodalIndex,
1223                                                ConstituentNodalValue);
1224
1225     delete[] ConstituentNodalIndex;
1226     delete[] ConstituentNodalValue;
1227
1228     delete[] ReverseDescendingConnectivityValue;
1229   }
1230   END_OF(LOC);
1231 }
1232
1233 //-----------------------------------------------------------------------------------------//
1234 //  void CONNECTIVITY::calculateReverseDescendingConnectivity(CONNECTIVITY *myConnectivity)
1235 //-----------------------------------------------------------------------------------------//
1236 //  {
1237 //    if (_entity==MED_CELL)
1238 //      throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
1239 //        Connectivity : we could not do that with MED_CELL entity !");
1240 //    if (myConnectivity->getEntity()!=MED_CELL)
1241 //      throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
1242 //              Connectivity : we need MED_CELL connectivity !");
1243 //    return;
1244 //  }
1245
1246 /*! Not implemented yet */
1247 //--------------------------------------------------------------------//
1248 void CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity)
1249 //--------------------------------------------------------------------//
1250 {
1251   // Mesh dimension !
1252   return;
1253 }
1254
1255
1256 /*!
1257   Give, for one element number of a specified entity the geometric type
1258   Of this element.
1259
1260   Example : medGeometryElement myType = myConnectivity.getElementType(MED_CELL,35);
1261 */
1262 //--------------------------------------------------------------------//
1263 medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity,int Number) const
1264 //--------------------------------------------------------------------//
1265 {
1266   if (_entity==Entity) {
1267     for(int i=1; i<=_numberOfTypes;i++)
1268       if (Number<_count[i])
1269         return _geometricTypes[i-1];
1270   }
1271   else if (_constituent!=NULL)
1272     return _constituent->getElementType(Entity,Number);
1273   else
1274     throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
1275   throw MEDEXCEPTION(LOCALIZED("getElementType : Wrong Number !"));
1276 }