src/MEDMEM/MEDMEM_Connectivity.cxx

   1 #include "MEDMEM_Connectivity.hxx"
   2 #include "MEDMEM_Family.hxx"
   3 #include "MEDMEM_Group.hxx"
   4 #include "MEDMEM_CellModel.hxx"
   5
   6 #include "MEDMEM_SkyLineArray.hxx"
   7 #include "MEDMEM_ModulusArray.hxx"
   8
   9 #include "MEDMEM_STRING.hxx"
  10 #include <iomanip>
  11
  12 using namespace std;
  13 using namespace MEDMEM;
  14 using namespace MED_EN;
  15
  16 // Enlarge the vector if necessary to insert the element
  17 static inline void insert_vector(vector<int> &Vect, int Indice, int Element)
  18 {
  19   if (Indice >= Vect.capacity())
  20     Vect.reserve(Indice + 1000);
  21
  22   if (Indice >= Vect.size())
  23     Vect.resize(Indice+1);
  24
  25   Vect[Indice] = Element;
  26 }
  27
  28 /*!
  29    Default Constructor. /n
  30    Default for Entity is MED_CELL and type of Connectivity is MED_NODAL */
  31 //--------------------------------------------------------------//
  32 CONNECTIVITY::CONNECTIVITY(medEntityMesh Entity /* =MED_CELL */) :
  33 //--------------------------------------------------------------//
  34                                 _entity(Entity),
  35                                 _typeConnectivity(MED_NODAL),
  36                                 _numberOfTypes(0),
  37                                 _geometricTypes((medGeometryElement*)NULL),
  38                                 _type((CELLMODEL*)NULL),
  39                                 _entityDimension(0),
  40                                 _numberOfNodes(0),
  41                                 _count((int*)NULL),
  42                                 _nodal((MEDSKYLINEARRAY*)NULL),
  43                                 _polygonsNodal((MEDSKYLINEARRAY*)NULL),
  44                                 _polyhedronNodal((POLYHEDRONARRAY*)NULL),
  45                                 _descending((MEDSKYLINEARRAY*)NULL),
  46                                 _polygonsDescending((MEDSKYLINEARRAY*)NULL),
  47                                 _polyhedronDescending((MEDSKYLINEARRAY*)NULL),
  48                                 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
  49                                 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
  50                                 _neighbourhood((MEDSKYLINEARRAY*)NULL),
  51                                 _constituent((CONNECTIVITY*)NULL)
  52 {
  53    MESSAGE("CONNECTIVITY(medEntityMesh Entity=MED_CELL)");
  54 }
  55
  56 /*!
  57    Constructor. /n
  58    Default for Entity is MED_CELL */
  59 //------------------------------------------------------------------------------//
  60 CONNECTIVITY::CONNECTIVITY(int numberOfTypes,medEntityMesh Entity /* =MED_CELL */):
  61 //------------------------------------------------------------------------------//
  62                                 _entity(Entity),
  63                                 _typeConnectivity(MED_NODAL),
  64                                 _numberOfTypes(numberOfTypes),
  65                                 _entityDimension(0),
  66                                 _numberOfNodes(0),
  67                                 _nodal((MEDSKYLINEARRAY*)NULL),
  68                                 _polygonsNodal((MEDSKYLINEARRAY*)NULL),
  69                                 _polyhedronNodal((POLYHEDRONARRAY*)NULL),
  70                                 _descending((MEDSKYLINEARRAY*)NULL),
  71                                 _polygonsDescending((MEDSKYLINEARRAY*)NULL),
  72                                 _polyhedronDescending((MEDSKYLINEARRAY*)NULL),
  73                                 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
  74                                 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
  75                                 _neighbourhood((MEDSKYLINEARRAY*)NULL),
  76                                 _constituent((CONNECTIVITY*)NULL)
  77 {
  78   MESSAGE("CONNECTIVITY(int numberOfTypes,medEntityMesh Entity=MED_CELL)");
  79   _geometricTypes = new medGeometryElement[numberOfTypes];
  80   _type = new CELLMODEL[numberOfTypes];
  81   _count = new int[numberOfTypes+1];
  82 }
  83
  84 /*!
  85   Copy Constructor.
  86 */
  87 //------------------------------------------------------//
  88 CONNECTIVITY::CONNECTIVITY (const CONNECTIVITY & m):
  89 //----------------------------------------------------//
  90                                 _entity           (m._entity),
  91                                 _typeConnectivity (m._typeConnectivity),
  92                                 _numberOfTypes    (m._numberOfTypes),
  93                                 _entityDimension  (m._entityDimension),
  94                                 _numberOfNodes    (m._numberOfNodes)
  95 {
  96  if (m._geometricTypes != NULL)
  97  {
  98     _geometricTypes = new medGeometryElement[_numberOfTypes];
  99     memcpy(_geometricTypes,m._geometricTypes,_numberOfTypes*sizeof(medGeometryElement));
 100  }
 101  else
 102     _geometricTypes = (medGeometryElement *) NULL;
 103
 104  if (m._type != NULL)
 105  {
 106     _type = new CELLMODEL[_numberOfTypes];
 107     for (int i=0;i<_numberOfTypes;i++)
 108         _type[i] = CELLMODEL(m._type[i]);
 109  }
 110  else
 111     _type = (CELLMODEL *) NULL;
 112
 113  if (m._count != NULL)
 114  {
 115       _count = new int[_numberOfTypes+1];
 116       memcpy(_count,m._count,(_numberOfTypes+1)*sizeof(int));
 117  }
 118  else
 119     _count = (int *) NULL;
 120
 121  if (m._nodal != NULL)
 122     _nodal = new MEDSKYLINEARRAY(* m._nodal);
 123  else
 124     _nodal = (MEDSKYLINEARRAY *) NULL;
 125
 126  if (m._polygonsNodal != NULL)
 127    _polygonsNodal = new MEDSKYLINEARRAY(* m._polygonsNodal);
 128  else
 129    _polygonsNodal = (MEDSKYLINEARRAY *) NULL;
 130
 131  if (m._polyhedronNodal != NULL)
 132    _polyhedronNodal = new POLYHEDRONARRAY(* m._polyhedronNodal);
 133  else
 134    _polyhedronNodal = (POLYHEDRONARRAY *) NULL;
 135
 136  if (m._descending != NULL)
 137     _descending = new MEDSKYLINEARRAY(* m._descending);
 138  else
 139     _descending = (MEDSKYLINEARRAY *) NULL;
 140
 141  if (m._polygonsDescending != NULL)
 142    _polygonsDescending = new MEDSKYLINEARRAY(* m._polygonsDescending);
 143  else
 144    _polygonsDescending = (MEDSKYLINEARRAY *) NULL;
 145
 146  if (m._polyhedronDescending != NULL)
 147    _polyhedronDescending = new MEDSKYLINEARRAY(* m._polyhedronDescending);
 148  else
 149    _polyhedronDescending = (MEDSKYLINEARRAY *) NULL;
 150
 151  if (m._reverseNodalConnectivity != NULL)
 152     _reverseNodalConnectivity = new MEDSKYLINEARRAY(* m._reverseNodalConnectivity);
 153  else
 154     _reverseNodalConnectivity = (MEDSKYLINEARRAY *) NULL;
 155
 156  if (m._reverseDescendingConnectivity != NULL)
 157     _reverseDescendingConnectivity = new MEDSKYLINEARRAY(* m._reverseDescendingConnectivity);
 158  else
 159     _reverseDescendingConnectivity = (MEDSKYLINEARRAY *) NULL;
 160
 161  if (m._neighbourhood != NULL)
 162     _neighbourhood = new MEDSKYLINEARRAY(* m._neighbourhood);
 163  else
 164     _neighbourhood = (MEDSKYLINEARRAY *) NULL;
 165
 166  if (m._constituent != NULL)
 167     _constituent = new CONNECTIVITY(* m._constituent);
 168  else
 169     _constituent = (CONNECTIVITY *) NULL;
 170 }
 171
 172 /*!
 173    Destructor./n
 174    desallocates existing pointers */
 175 //----------------------------//
 176 CONNECTIVITY::~CONNECTIVITY()
 177 //----------------------------//
 178 {
 179   MESSAGE("Destructeur de CONNECTIVITY()");
 180
 181   if (_geometricTypes != NULL)
 182      delete [] _geometricTypes;
 183   if (_type != NULL)
 184      delete [] _type;
 185   if (_count != NULL)
 186      delete [] _count;
 187   if (_nodal != NULL)
 188      delete _nodal;
 189   if (_polygonsNodal != NULL)
 190      delete _polygonsNodal;
 191   if (_polyhedronNodal != NULL)
 192      delete _polyhedronNodal;
 193   if (_descending != NULL)
 194      delete _descending;
 195   if (_polygonsDescending != NULL)
 196      delete _polygonsDescending;
 197   if (_polyhedronDescending != NULL)
 198     delete _polyhedronDescending;
 199   if (_reverseNodalConnectivity != NULL)
 200      delete _reverseNodalConnectivity;
 201   if (_reverseDescendingConnectivity != NULL)
 202      delete _reverseDescendingConnectivity;
 203   if (_constituent != NULL)
 204      delete _constituent;
 205 }
 206
 207 /*!
 208    set _constituent to Constituent
 209    be aware desallocation of _constituent is done by CONNECTIVITY:~CONNECTIVITY
 210    throws an exception if Constituent = MED_CELL
 211    A DOCUMENTER
 212     */
 213 //----------------------------------------------------------//
 214 void CONNECTIVITY::setConstituent(CONNECTIVITY * Constituent)
 215                     throw (MEDEXCEPTION)
 216 //----------------------------------------------------------//
 217 {
 218   medEntityMesh Entity = Constituent->getEntity();
 219   if (Entity == MED_CELL)
 220     throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : could not set constituent on MED_CELL !"));
 221
 222   if ((Entity == MED_EDGE)&(_entityDimension == 3))
 223   {
 224     if (_constituent == NULL)
 225       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : Entity not found !"));
 226     _constituent->setConstituent(Constituent);
 227   }
 228   else
 229     _constituent = Constituent;
 230 }
 231
 232 /*! Duplicated Types array in CONNECTIVITY object. */
 233 //--------------------------------------------------------------------//
 234 void CONNECTIVITY::setGeometricTypes(const medGeometryElement * Types,
 235                                      const medEntityMesh Entity)
 236                                      throw (MEDEXCEPTION)
 237 //--------------------------------------------------------------------//
 238 {
 239   if (Entity == _entity)
 240     for (int i=0; i<_numberOfTypes; i++)
 241     {
 242       _geometricTypes[i] = Types[i];
 243       _type[i] = CELLMODEL(_geometricTypes[i]);
 244     }
 245   else
 246     {
 247     if (_constituent == NULL)
 248         throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setGeometricTypes : Entity not found !"));
 249     _constituent->setGeometricTypes(Types,Entity);
 250   }
 251 }
 252
 253 /*! A DOCUMENTER */
 254 //--------------------------------------------------------------------//
 255 void CONNECTIVITY::setCount(const int * Count, const medEntityMesh Entity)
 256                             throw (MEDEXCEPTION)
 257 //--------------------------------------------------------------------//
 258 {
 259   if (Entity == _entity)
 260   {
 261     int * index = new int[Count[_numberOfTypes]];
 262     index[0]=1;
 263     _count[0]=1;
 264     for (int i=0; i<_numberOfTypes; i++) {
 265       _count[i+1] = Count[i+1];
 266       int NumberOfNodesPerElement = _type[i].getNumberOfNodes() ;
 267       for (int j=_count[i]; j<_count[i+1]; j++)
 268         index[j] = index[j-1]+NumberOfNodesPerElement;
 269     }
 270     // allocate _nodal
 271     if (_nodal != NULL) delete _nodal;
 272     _nodal = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,index[_count[_numberOfTypes]-1]-1);
 273     _nodal->setIndex(index);
 274     delete [] index;
 275   }
 276   else
 277   {
 278     if (_constituent == NULL)
 279       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setCount : Entity not found !"));
 280     _constituent->setCount(Count,Entity);
 281   }
 282 }
 283
 284 //--------------------------------------------------------------------//
 285 void CONNECTIVITY::setNodal(const int * Connectivity,
 286                             const medEntityMesh Entity,
 287                             const medGeometryElement Type)
 288                             throw (MEDEXCEPTION)
 289 //--------------------------------------------------------------------//
 290 {
 291   if (_entity == Entity)
 292   {
 293     // find geometric type
 294     bool find = false;
 295     for (int i=0; i<_numberOfTypes; i++)
 296       if (_geometricTypes[i] == Type) {
 297         find = true;
 298         int NumberOfNodePerElements = _type[i].getNumberOfNodes() ;
 299         //_nodal->setI(i+1,Connectivity);
 300         for( int j=_count[i];j<_count[i+1]; j++)
 301           _nodal->setI(j,Connectivity+(j-_count[i])*NumberOfNodePerElements);
 302       }
 303     if (!find)
 304       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : geometric type not found"));
 305   } else
 306   {
 307     if (_constituent == NULL)
 308       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : Entity not found !"));
 309       _constituent->setNodal(Connectivity,Entity,Type);
 310   }
 311 }
 312
 313
 314 //--------------------------------------------------------------------//
 315 void CONNECTIVITY::setPolygonsConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, const int* PolygonsConnectivity, const int* PolygonsConnectivityIndex, int ConnectivitySize, int NumberOfPolygons)
 316 //--------------------------------------------------------------------//
 317 {
 318   const char* LOC = "CONNECTIVITY::setPolygonsConnectivity";
 319   BEGIN_OF(LOC);
 320
 321   if (_entity == Entity)
 322     {
 323       MEDSKYLINEARRAY* Connectivity = new MEDSKYLINEARRAY(NumberOfPolygons,ConnectivitySize,PolygonsConnectivityIndex,PolygonsConnectivity);
 324
 325       if (ConnectivityType == MED_NODAL)
 326         {
 327           if (_polygonsNodal != (MEDSKYLINEARRAY*) NULL)
 328             delete _polygonsNodal;
 329           _polygonsNodal = Connectivity;
 330         }
 331       else
 332         {
 333           if (_typeConnectivity != MED_DESCENDING)
 334             _typeConnectivity = MED_DESCENDING; //by default it is set to MED_NODAL
 335           if (_polygonsDescending != (MEDSKYLINEARRAY*) NULL)
 336             delete _polygonsDescending;
 337           _polygonsDescending = Connectivity;
 338         }
 339     }
 340   else
 341     {
 342       if (_constituent == (CONNECTIVITY*) NULL)
 343         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not found !"));
 344       _constituent->setPolygonsConnectivity(ConnectivityType, Entity, PolygonsConnectivity, PolygonsConnectivityIndex, ConnectivitySize, NumberOfPolygons);
 345     }
 346 }
 347
 348
 349 //--------------------------------------------------------------------//
 350 void CONNECTIVITY::setPolyhedronConnectivity(medConnectivity ConnectivityType, const int* PolyhedronConnectivity, const int* PolyhedronIndex, int ConnectivitySize, int NumberOfPolyhedron, const int* PolyhedronFacesIndex /* =(int*)NULL */, int NumberOfFaces /* =0 */)
 351 //--------------------------------------------------------------------//
 352 {
 353   const char* LOC = "CONNECTIVITY::setPolyhedronConnectivity";
 354   BEGIN_OF(LOC);
 355
 356   if (_entity == MED_CELL)
 357     {
 358       if (ConnectivityType == MED_NODAL)
 359         {
 360           if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
 361             delete _polyhedronNodal;
 362           _polyhedronNodal = new POLYHEDRONARRAY(NumberOfPolyhedron,NumberOfFaces,ConnectivitySize);
 363           _polyhedronNodal->setPolyhedronIndex(PolyhedronIndex);
 364           _polyhedronNodal->setFacesIndex(PolyhedronFacesIndex);
 365           _polyhedronNodal->setNodes(PolyhedronConnectivity);
 366         }
 367       else
 368         {
 369           if (_typeConnectivity != MED_DESCENDING)
 370             _typeConnectivity = MED_DESCENDING; //by default it is set to MED_NODAL
 371           if (_polyhedronDescending != (MEDSKYLINEARRAY*) NULL)
 372             delete _polyhedronDescending;
 373           _polyhedronDescending = new MEDSKYLINEARRAY(NumberOfPolyhedron,ConnectivitySize,PolyhedronIndex,PolyhedronConnectivity);
 374         }
 375     }
 376   else
 377     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : _entity must be MED_CELL to set polyhedron !"));
 378 }
 379
 380
 381 /*! A DOCUMENTER */
 382 //------------------------------------------------------------------------------------------//
 383 void CONNECTIVITY::calculateConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity)
 384 //------------------------------------------------------------------------------------------//
 385 {
 386   MESSAGE("CONNECTIVITY::calculateConnectivity");
 387
 388   // a temporary limitation due to calculteDescendingConnectivity function !!!
 389   if ((_entityDimension==3) & (Entity==MED_EDGE))
 390     throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build EDGE !");
 391
 392   if (Entity==_entity)
 393     if (ConnectivityType==MED_NODAL)
 394       calculateNodalConnectivity();
 395     else
 396       if (Entity==MED_CELL)
 397         calculateDescendingConnectivity();
 398       else
 399         throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build DESCENDING !");
 400   if (Entity!=_entity) {
 401     calculateDescendingConnectivity();
 402     if (_entityDimension == 2 || _entityDimension == 3)
 403       _constituent->calculateConnectivity(ConnectivityType,Entity);
 404   }
 405 }
 406
 407 /*!  Give, in full or no interlace mode (for nodal connectivity),
 408      descending or nodal connectivity.
 409
 410       You must give a %medEntityMesh (ie:MED_EDGE)
 411       and a %medGeometryElement (ie:MED_SEG3).  */
 412
 413 //------------------------------------------------------------//
 414 void CONNECTIVITY::updateFamily(const vector<FAMILY*>& myFamilies)
 415 //------------------------------------------------------------//
 416 {
 417   const char * LOC = "CONNECTIVITY::updateFamily(vector<FAMILY*>) ";
 418   int numberOfFamilies = myFamilies.size();
 419   if (numberOfFamilies == 0 || _constituent == NULL)
 420     return;
 421   // does we do an update ?
 422   if ((_constituent != NULL) && (_descending != NULL))
 423     return;
 424   if (myFamilies[0]->getEntity() != _constituent->getEntity())
 425     return;
 426   CONNECTIVITY * oldConstituent = _constituent;
 427   _constituent = (CONNECTIVITY *)NULL;
 428   if (oldConstituent->_nodal==NULL)
 429     throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
 430
 431   //Loc vars defined to treat polygons exactly the same as classic types. Not nice but necessary.
 432   int oldNumberOfFaceTab[2];
 433   const int * oldConstituentValueTab[2];
 434   const int * oldConstituentIndexTab[2];
 435   int * renumberingFromOldToNewTab[2];//Final mapping array between old numbers and new numbers.;
 436
 437   int oldNumberOfFace = oldConstituent->_nodal->getNumberOf(); oldNumberOfFaceTab[0]=oldNumberOfFace;
 438   const int * oldConstituentValue = oldConstituent->_nodal->getValue(); oldConstituentValueTab[0]=oldConstituentValue;
 439   const int * oldConstituentIndex = oldConstituent->_nodal->getIndex(); oldConstituentIndexTab[0]=oldConstituentIndex;
 440   int * renumberingFromOldToNew= new int [oldNumberOfFace]; renumberingFromOldToNewTab[0]=renumberingFromOldToNew;
 441
 442   int oldNumberOfFacePoly = oldConstituent->getNumberOfPolygons();
 443   const int * oldConstituentValuePoly=0;
 444   const int * oldConstituentIndexPoly=0;
 445   int * renumberingFromOldToNewPoly=0;
 446
 447   int nbOfTurnInGlobalLoop=1;//Defined to know if a second search on polygons is needed.
 448   if(oldNumberOfFacePoly>0)
 449     {
 450       oldNumberOfFaceTab[1]=oldNumberOfFacePoly;
 451       oldConstituentValuePoly = oldConstituent->_polygonsNodal->getValue(); oldConstituentValueTab[1]=oldConstituentValuePoly;
 452       oldConstituentIndexPoly = oldConstituent->_polygonsNodal->getIndex(); oldConstituentIndexTab[1]=oldConstituentIndexPoly;
 453       renumberingFromOldToNewPoly=new int[oldNumberOfFacePoly]; renumberingFromOldToNewTab[1]=renumberingFromOldToNewPoly;
 454       nbOfTurnInGlobalLoop++;
 455     }
 456
 457   calculateDescendingConnectivity();//perform calculation of descending connectivity to update all connectivities after having taken into account of direction of faces...
 458   _constituent->calculateReverseNodalConnectivity();//getting d-1 nodal connectivity to find get new face numbers from nodes numbers...
 459
 460   const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity(); //Common to polygons and classic geometric types
 461   const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex(); //Common to polygons and classic geometric types
 462
 463   for(int loop=0;loop<nbOfTurnInGlobalLoop;loop++)
 464     {
 465       int oldNumberOfFaceLoop=oldNumberOfFaceTab[loop];
 466       const int * oldConstituentValueLoop=oldConstituentValueTab[loop];
 467       const int * oldConstituentIndexLoop= oldConstituentIndexTab[loop];
 468       int * renumberingFromOldToNewLoop=renumberingFromOldToNewTab[loop];
 469       for(int iOldFace=0;iOldFace<oldNumberOfFaceLoop;iOldFace++)
 470         {
 471           const int *nodesOfCurrentFaceOld=oldConstituentValueLoop+oldConstituentIndexLoop[iOldFace]-1;
 472           int nbOfNodesOfCurrentFaceOld=oldConstituentIndexLoop[iOldFace+1]-oldConstituentIndexLoop[iOldFace];
 473
 474           //retrieving new number of polygon face...
 475           int sizeOfNewFaceNb=reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]]-reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]-1];
 476           int *newFaceNb=new int[sizeOfNewFaceNb];
 477           memcpy(newFaceNb,reverseFaceNodal+reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]-1]-1,sizeOfNewFaceNb*sizeof(int));//all faces candidates are in newFaceNb !!!
 478           for(int curNode=1;curNode<nbOfNodesOfCurrentFaceOld && sizeOfNewFaceNb>1;curNode++)
 479             {
 480               const int *newFacesNbForCurNode=reverseFaceNodal+reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]-1]-1;
 481               int sizeOfNewFacesNbForCurNode=reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]]-reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]-1];
 482               for(int i=0;i<sizeOfNewFaceNb;)
 483                 {
 484                   bool found=false;
 485                   for(int j=0;j<sizeOfNewFacesNbForCurNode && !found;j++)
 486                     {
 487                       if(newFacesNbForCurNode[j]==newFaceNb[i])
 488                         found=true;
 489                     }
 490                   if(found)
 491                     i++;
 492                   else
 493                     newFaceNb[i]=newFaceNb[--sizeOfNewFaceNb];
 494                 }
 495             }
 496           if(sizeOfNewFaceNb!=1)
 497             throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"unexisting face specified!!!"));
 498           renumberingFromOldToNewLoop[iOldFace]=newFaceNb[0];
 499           delete [] newFaceNb;
 500           //end of retrieving new number of polygon face...
 501           int nbOfNodesOfCurrentFaceNew;
 502           const int *nodesOfCurrentFaceNew=_constituent->getConnectivityOfAnElementWithPoly(MED_NODAL,_constituent->getEntity(),
 503                                                                                             renumberingFromOldToNewLoop[iOldFace],nbOfNodesOfCurrentFaceNew);
 504           MEDMODULUSARRAY modulusArrayOld(nbOfNodesOfCurrentFaceOld,nodesOfCurrentFaceOld);
 505           MEDMODULUSARRAY modulusArrayNew(nbOfNodesOfCurrentFaceNew,nodesOfCurrentFaceNew);
 506           int retCompareNewOld=modulusArrayNew.compare(modulusArrayOld);
 507           if(retCompareNewOld==0)
 508             throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"Uncompatible given user face with calculated existing faces"));
 509           if(retCompareNewOld==-1)
 510             invertConnectivityForAFace(renumberingFromOldToNewLoop[iOldFace],nodesOfCurrentFaceOld,loop==1);
 511         }
 512     }
 513   // Updating the Family
 514   for(vector<FAMILY*>::const_iterator iter=myFamilies.begin();iter!=myFamilies.end();iter++)
 515     (*iter)->changeElementsNbs(_constituent->getEntity(),renumberingFromOldToNew,oldNumberOfFace,renumberingFromOldToNewPoly);
 516   delete oldConstituent ;
 517   delete [] renumberingFromOldToNew;
 518   if(oldNumberOfFacePoly>0)
 519     delete [] renumberingFromOldToNewPoly;
 520   return;
 521 }
 522
 523 //------------------------------------------------------------------------------------------------------------------//
 524 const int * MEDMEM::CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type)
 525 //------------------------------------------------------------------------------------------------------------------//
 526 {
 527   const char * LOC = "CONNECTIVITY::getConnectivity";
 528   //  BEGIN_OF(LOC);
 529
 530   MEDSKYLINEARRAY * Connectivity;
 531   if (Entity==_entity) {
 532
 533     if (ConnectivityType==MED_NODAL)
 534       {
 535         calculateNodalConnectivity();
 536         Connectivity=_nodal;
 537       }
 538     else
 539       {
 540         calculateDescendingConnectivity();
 541         Connectivity=_descending;
 542       }
 543
 544     if (Connectivity!=NULL)
 545       if (Type==MED_ALL_ELEMENTS)
 546         return Connectivity->getValue();
 547       else {
 548         for (int i=0; i<_numberOfTypes; i++)
 549           if (_geometricTypes[i]==Type)
 550             //return Connectivity->getI(i+1);
 551             return Connectivity->getI(_count[i]);
 552         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
 553       }
 554     else
 555       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
 556   } else
 557     if (_constituent != NULL)
 558       return _constituent->getConnectivity(ConnectivityType,Entity,Type);
 559
 560   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 561 }
 562
 563 //------------------------------------------------------------------------------------------------------------------//
 564 int CONNECTIVITY::getConnectivityLength(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type)
 565 //------------------------------------------------------------------------------------------------------------------//
 566 {
 567   const char * LOC = "CONNECTIVITY::getConnectivity";
 568   BEGIN_OF(LOC);
 569
 570   MEDSKYLINEARRAY * Connectivity;
 571   if (Entity==_entity) {
 572
 573     if (ConnectivityType==MED_NODAL)
 574       {
 575         calculateNodalConnectivity();
 576         Connectivity=_nodal;
 577       }
 578     else
 579       {
 580         calculateDescendingConnectivity();
 581         Connectivity=_descending;
 582       }
 583
 584     if (Connectivity!=NULL)
 585       if (Type==MED_ALL_ELEMENTS)
 586         return Connectivity->getLength();
 587       else {
 588         for (int i=0; i<_numberOfTypes; i++)
 589           if (_geometricTypes[i]==Type)
 590             return Connectivity->getNumberOfI(_count[i]);
 591         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
 592       }
 593     else
 594       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
 595   }
 596   else
 597     if (_constituent != NULL)
 598       return _constituent->getConnectivityLength(ConnectivityType,Entity,Type);
 599
 600   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 601 }
 602
 603 /*!  Give morse index array to use with
 604      getConnectivity(MED_FULL_INTERLACE,mode,entity,MED_ALL_ELEMENTS).
 605
 606       Each value give start index for corresponding entity in connectivity array.
 607
 608       Example : i-th element, j-th node of it :
 609       - In C mode : Connectivity[ConnectivityIndex[i]-1+j-1]
 610       - In fortran mode : Connectivity[ConnectivityIndex[i]+j] */
 611 //-----------------------------------------------------------------------------------------------//
 612 const int * CONNECTIVITY::getConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity)
 613 //-----------------------------------------------------------------------------------------------//
 614 {
 615   const char * LOC = "CONNECTIVITY::getConnectivityIndex";
 616   BEGIN_OF(LOC);
 617
 618   MEDSKYLINEARRAY * Connectivity;
 619   if (Entity==_entity) {
 620
 621     if (ConnectivityType==MED_NODAL)
 622       Connectivity=_nodal;
 623     else
 624       Connectivity=_descending;
 625
 626     if (Connectivity!=NULL)
 627       return Connectivity->getIndex();
 628     else
 629       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
 630
 631   } else
 632     if (_constituent != NULL)
 633       return _constituent->getConnectivityIndex(ConnectivityType,Entity);
 634
 635   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 636 }
 637
 638
 639 //-------------------------------------------------------------//
 640 const int* CONNECTIVITY::getPolygonsConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity)
 641 //-------------------------------------------------------------//
 642 {
 643   const char* LOC = "CONNECTIVITY::getPolygonsConnectivity";
 644   BEGIN_OF(LOC);
 645
 646   MEDSKYLINEARRAY* Connectivity;
 647   if (Entity == _entity)
 648     {
 649       if (ConnectivityType == MED_NODAL)
 650         {
 651           calculateNodalConnectivity();
 652           Connectivity = _polygonsNodal;
 653         }
 654       else
 655         {
 656           calculateDescendingConnectivity();
 657           Connectivity = _polygonsDescending;
 658         }
 659       if (Connectivity != NULL)
 660         return Connectivity->getValue();
 661       else
 662         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polygons Connectivity not defined !"));
 663     }
 664   else
 665     if (_constituent != NULL)
 666       return _constituent->getPolygonsConnectivity(ConnectivityType, Entity);
 667   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 668 }
 669
 670
 671 //-------------------------------------------------------------//
 672 const int* CONNECTIVITY::getPolygonsConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity)
 673 //-------------------------------------------------------------//
 674 {
 675   const char* LOC = "CONNECTIVITY::getPolygonsConnectivityIndex";
 676   BEGIN_OF(LOC);
 677
 678   MEDSKYLINEARRAY* Connectivity;
 679   if (Entity == _entity)
 680     {
 681       if (ConnectivityType == MED_NODAL)
 682         {
 683           //      calculateNodalConnectivity();
 684           Connectivity = _polygonsNodal;
 685         }
 686       else
 687         {
 688           //      calculateDescendingConnectivity();
 689           Connectivity = _polygonsDescending;
 690         }
 691       if (Connectivity != NULL)
 692         return Connectivity->getIndex();
 693       else
 694         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polygons Connectivity not defined !"));
 695     }
 696   else
 697     if (_constituent != NULL)
 698       return _constituent->getPolygonsConnectivityIndex(ConnectivityType, Entity);
 699   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 700 }
 701
 702
 703 /*! We suppose in this method that nodal and descending connectivities
 704   are coherent.*/
 705 //-------------------------------------------------------------//
 706 int CONNECTIVITY::getNumberOfPolygons() const
 707 //-------------------------------------------------------------//
 708 {
 709   if (_polygonsNodal != (MEDSKYLINEARRAY*) NULL)
 710     return _polygonsNodal->getNumberOf();
 711   else if (_polygonsDescending != (MEDSKYLINEARRAY*) NULL)
 712     return _polygonsDescending->getNumberOf();
 713   else
 714     return 0;
 715 }
 716
 717
 718 //--------------------------------------------------------------//
 719 const int* CONNECTIVITY::getPolyhedronConnectivity(medConnectivity ConnectivityType) const
 720 //--------------------------------------------------------------//
 721 {
 722   const char* LOC = "CONNECTIVITY::getPolyhedronConnectivity";
 723   BEGIN_OF(LOC);
 724
 725   if (_entity == MED_CELL) //polyhedron can only be MED_CELL
 726     {
 727       if (ConnectivityType == MED_NODAL)
 728         {
 729           ((CONNECTIVITY *)(this))->calculateNodalConnectivity();
 730           if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
 731             return _polyhedronNodal->getNodes();
 732           else
 733             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polyhedron Nodal Connectivity not defined !"));
 734         }
 735       else
 736         {
 737           ((CONNECTIVITY *)(this))->calculateDescendingConnectivity();
 738           if (_polyhedronDescending != (MEDSKYLINEARRAY*) NULL)
 739             return _polyhedronDescending->getValue();
 740           else
 741             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polyhedron Descending Connectivity not defined !"));
 742         }
 743     }
 744   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : No polyhedron for _entity different from MED_CELL !"));
 745 }
 746
 747
 748 //---------------------------------------------------------------//
 749 const int* CONNECTIVITY::getPolyhedronFacesIndex() const
 750 //---------------------------------------------------------------//
 751 {
 752   const char* LOC = "CONNECTIVITY::getPolyhedronFacesIndex";
 753   BEGIN_OF(LOC);
 754
 755   if (_entity == MED_CELL) //polyhedron can only be MED_CELL
 756     {
 757       //      calculateNodalConnectivity();
 758       if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
 759         return _polyhedronNodal->getFacesIndex();
 760       else
 761         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : No Polyhedron in that Connectivity !"));
 762     }
 763   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : No polyhedron for _entity different from MED_CELL !"));
 764 }
 765
 766
 767 //--------------------------------------------------------------//
 768 const int* CONNECTIVITY::getPolyhedronIndex(medConnectivity ConnectivityType) const
 769 //--------------------------------------------------------------//
 770 {
 771   const char* LOC = "CONNECTIVITY::getPolyhedronIndex";
 772   BEGIN_OF(LOC);
 773
 774   if (_entity == MED_CELL) //polyhedron can only be MED_CELL
 775     {
 776       if (ConnectivityType == MED_NODAL)
 777         {
 778           //      calculateNodalConnectivity();
 779           if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
 780             return _polyhedronNodal->getPolyhedronIndex();
 781           else
 782             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polyhedron Nodal Connectivity not defined !"));
 783         }
 784       else
 785         {
 786           //      calculateDescendingConnectivity();
 787           if (_polyhedronDescending != (MEDSKYLINEARRAY*) NULL)
 788             return _polyhedronDescending->getIndex();
 789           else
 790             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polyhedron Descending Connectivity not defined !"));
 791         }
 792     }
 793   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : No polyhedron for _entity different from MED_CELL !"));
 794 }
 795
 796
 797 /*! We suppose in this method that nodal and descending connectivities
 798   are coherent.*/
 799 //-------------------------------------------------------------//
 800 int CONNECTIVITY::getNumberOfPolyhedronFaces() const
 801 //-------------------------------------------------------------//
 802 {
 803   //  if (_polyhedronNodal == (POLYHEDRONARRAY*) NULL)
 804   //    calculateNodalConnectivity();
 805   if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
 806     return _polyhedronNodal->getNumberOfFaces();
 807   else
 808     return 0;
 809 }
 810
 811
 812 /*! We suppose in this method that nodal and descending connectivities
 813   are coherent.*/
 814 //--------------------------------------------------------------//
 815 int CONNECTIVITY::getNumberOfPolyhedron() const
 816 //--------------------------------------------------------------//
 817 {
 818   if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
 819     return _polyhedronNodal->getNumberOfPolyhedron();
 820   else if (_polyhedronDescending != (MEDSKYLINEARRAY*) NULL)
 821     return _polyhedronDescending->getNumberOf();
 822   else
 823     return 0;
 824 }
 825
 826
 827 /*! A DOCUMENTER */
 828 //--------------------------------------------------------------//
 829 const CELLMODEL & CONNECTIVITY::getType(medGeometryElement Type) const
 830 //--------------------------------------------------------------//
 831 {
 832   const char * LOC = "CONNECTIVITY::getType";
 833   BEGIN_OF(LOC);
 834
 835   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
 836     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE !"));
 837   for (int i=0; i<_numberOfTypes; i++)
 838     if (_geometricTypes[i]==Type)
 839       return _type[i];
 840   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
 841 }
 842
 843 /*!  Returns the number of elements of type %medGeometryElement.
 844      Note : not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
 845 //------------------------------------------------------------------------//
 846 int CONNECTIVITY::getNumberOfNodesInType(medGeometryElement Type) const
 847 //------------------------------------------------------------------------//
 848 {
 849   const char * LOC = "CONNECTIVITY::getNumberOfNodesInType";
 850   BEGIN_OF(LOC);
 851
 852   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
 853     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!"));
 854   for (int i=0; i<_numberOfTypes; i++)
 855     if (_geometricTypes[i]==Type)
 856       return _type[i].getNumberOfNodes();
 857   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
 858 }
 859
 860 /*!  Returns the number of geometric sub cells of %medGeometryElement type.
 861 not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
 862 //------------------------------------------------------------------------//
 863 int CONNECTIVITY::getNumberOfSubCellInType(medGeometryElement Type) const
 864 //------------------------------------------------------------------------//
 865 {
 866   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
 867     throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!");
 868   for (int i=0; i<_numberOfTypes; i++)
 869     if (_geometricTypes[i]==Type)
 870       return _type[i].getNumberOfConstituents(1);
 871   throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement not found !");
 872 }
 873
 874 /*!  Returns the number of elements of type %medGeometryElement.
 875
 876      Note :
 877       - Implemented for MED_ALL_ELEMENTS
 878       - Not implemented for MED_ALL_ENTITIES (A VERIFIER)
 879       - Not implemented for MED_NONE */
 880 //-----------------------------------------------------------------------------------//
 881 int CONNECTIVITY::getNumberOf(medEntityMesh Entity, medGeometryElement Type) const
 882 //-----------------------------------------------------------------------------------//
 883 {
 884   //const char * LOC = "CONNECTIVITY::getNumberOf";
 885   if (Entity==_entity) {
 886     if (Type==MED_NONE)
 887       return 0; // not defined !
 888     //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_NONE"));
 889     if (!existConnectivity(MED_NODAL,Entity) && !existConnectivity(MED_DESCENDING,Entity)) return 0; //case with only polygons for example
 890     if (Type==MED_ALL_ELEMENTS)
 891       return _count[_numberOfTypes]-1;
 892     for (int i=0; i<_numberOfTypes; i++)
 893       if (_geometricTypes[i]==Type)
 894         return (_count[i+1] - _count[i]);
 895     //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
 896   } else
 897     if (_constituent != NULL)
 898       return _constituent->getNumberOf(Entity,Type);
 899
 900   return 0; // valid if they are nothing else !
 901   //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 902 }
 903
 904 /*! A DOCUMENTER */
 905 //--------------------------------------------------------------//
 906 const int* CONNECTIVITY::getValue(medConnectivity TypeConnectivity,
 907                                 medGeometryElement Type)
 908 //--------------------------------------------------------------//
 909 {
 910   if (TypeConnectivity == MED_NODAL)
 911     {
 912           calculateNodalConnectivity();
 913           if (Type==MED_ALL_ELEMENTS)
 914             return _nodal->getValue();
 915           for (int i=0; i<_numberOfTypes; i++)
 916             if (_geometricTypes[i]==Type)
 917               return _nodal->getI(_count[i]);
 918     }
 919   else
 920     {
 921       calculateDescendingConnectivity();
 922       if (Type==MED_ALL_ELEMENTS)
 923         return _descending->getValue();
 924       for (int i=0; i<_numberOfTypes; i++)
 925         if (_geometricTypes[i]==Type)
 926           return _descending->getI(Type);
 927     }
 928   throw MEDEXCEPTION("Not found");
 929 }
 930
 931 /*! A DOCUMENTER */
 932 //---------------------------------------------------------------------//
 933 const int* CONNECTIVITY:: getValueIndex(medConnectivity TypeConnectivity)
 934 //---------------------------------------------------------------------//
 935 {
 936   if (TypeConnectivity == MED_NODAL)
 937     {
 938       calculateNodalConnectivity();
 939       return _nodal->getIndex();
 940     }
 941   else
 942     {
 943       calculateDescendingConnectivity();
 944       return _descending->getIndex();
 945     }
 946 }
 947
 948 /*! Not yet implemented */
 949 //----------------------------------------------//
 950 const int* CONNECTIVITY:: getNeighbourhood() const
 951 //----------------------------------------------//
 952 {
 953   throw MEDEXCEPTION("CONNECTIVITY::getNeighbourhood : Not implemented");
 954 }
 955
 956 /*! Returns an array which contains, for each node, all cells
 957     arround it. */
 958 //-------------------------------------------------//
 959 const int* CONNECTIVITY::getReverseNodalConnectivity()
 960 //-------------------------------------------------//
 961 {
 962   calculateReverseNodalConnectivity();
 963   return _reverseNodalConnectivity->getValue();
 964 }
 965
 966 /*!  Give index array to use with getReverseConnectivity(MED_NODAL).
 967      It is unusefull with MED_DESCENDING mode, because we have allways two cells.  */
 968 //-------------------------------------------------------//
 969 const int* CONNECTIVITY::getReverseNodalConnectivityIndex()
 970 //-------------------------------------------------------//
 971 {
 972   calculateReverseNodalConnectivity();
 973   return _reverseNodalConnectivity->getIndex();
 974 }
 975
 976 /*! Returns an array which contains, for each face (or edge),
 977     the 2 cells of each side. First is cell which face normal is outgoing.
 978     arround it. */
 979 //------------------------------------------------------//
 980 const int* CONNECTIVITY::getReverseDescendingConnectivity()
 981 //------------------------------------------------------//
 982 {
 983   // it is in _constituent connectivity only if we are in MED_CELL
 984   // (we could not for instance calculate face-edge connectivity !)
 985   if (_entity!=MED_CELL)
 986     throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivity : Error Only in MED_CELL connectivity");
 987
 988   // we need descending connectivity
 989   calculateDescendingConnectivity();
 990   return _reverseDescendingConnectivity->getValue();
 991 }
 992
 993 /*! calculate the reverse descending Connectivity
 994     and returns the index  ( A DOCUMENTER MIEUX)*/
 995 //-----------------------------------------------------------//
 996 const int* CONNECTIVITY::getReverseDescendingConnectivityIndex()
 997 //-----------------------------------------------------------//
 998 {
 999   // it is in _constituent connectivity only if we are in MED_CELL
1000   if (_entity!=MED_CELL)
1001     throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivityIndex : Error Only in MED_CELL connectivity");
1002
1003   // we need descending connectivity
1004   calculateDescendingConnectivity();
1005   return _reverseDescendingConnectivity->getIndex();
1006 }
1007
1008 /*! A DOCUMENTER (et a finir ???) */
1009 //--------------------------------------------//
1010 void CONNECTIVITY::calculateNodalConnectivity()
1011 //--------------------------------------------//
1012 {
1013   if (_nodal==NULL && _polygonsNodal==NULL && _polyhedronNodal==NULL)
1014     {
1015       if (_descending==NULL && _polygonsDescending==NULL && _polyhedronDescending==NULL)
1016         throw MEDEXCEPTION("CONNECTIVITY::calculateNodalConnectivity : No connectivity found !");
1017       // calculate _nodal, _polygonsNodal and _polyhedronNodal from _descending, _polygonsDescending and _polyhedronDescending
1018     }
1019 }
1020
1021 /*! If not yet done, calculate the nodal Connectivity
1022     and the reverse nodal Connectivity*/
1023 //---------------------------------------------------//
1024 void CONNECTIVITY::calculateReverseNodalConnectivity()
1025 //---------------------------------------------------//
1026 {
1027   const char * LOC = "CONNECTIVITY::calculateReverseNodalConnectivity : ";
1028   BEGIN_OF(LOC);
1029
1030   SCRUTE(_nodal);
1031   SCRUTE(_reverseNodalConnectivity);
1032
1033   if (_nodal==NULL)
1034     calculateNodalConnectivity();
1035
1036   if(_reverseNodalConnectivity==NULL) {
1037
1038     int node_number = 0;
1039     vector <vector <int> > reverse_connectivity;
1040     reverse_connectivity.resize(_numberOfNodes+1);
1041
1042   // Treat all cells types
1043
1044     for (int j = 0; j < _numberOfTypes; j++)
1045       {
1046         // node number of the cell type
1047         node_number = _type[j].getNumberOfNodes();
1048         // treat all cells of a particular type
1049         for (int k = _count[j]; k < _count[j+1]; k++)
1050           // treat all nodes of the cell type
1051           for (int local_node_number = 1; local_node_number < node_number+1; local_node_number++)
1052             {
1053               int global_node = _nodal->getIJ(k,local_node_number);
1054               reverse_connectivity[global_node].push_back(k);
1055             }
1056       }
1057
1058     if(_polygonsNodal != NULL && _entityDimension==2)
1059       {
1060         int nbOfPolygons=_polygonsNodal->getNumberOf();
1061         int offset=getNumberOf(_entity,MED_ALL_ELEMENTS);
1062         const int *index=_polygonsNodal->getIndex();
1063         const int *value=_polygonsNodal->getValue();
1064         for (int local_polyg_number = 0; local_polyg_number < nbOfPolygons; local_polyg_number++)
1065           {
1066             for(int i=index[local_polyg_number];i<index[local_polyg_number+1];i++)
1067               reverse_connectivity[value[i-1]].push_back(offset+local_polyg_number+1);
1068           }
1069       }
1070
1071     if(_polyhedronNodal != NULL && _entityDimension==3)
1072       {
1073         int nbOfPolyhedra=_polyhedronNodal->getNumberOfPolyhedron();
1074         int offset=getNumberOf(_entity,MED_ALL_ELEMENTS);
1075         for (int local_polyh_number = 0; local_polyh_number < nbOfPolyhedra; local_polyh_number++)
1076             {
1077               int nbOfNodes;
1078               int global_polyh_number=offset+local_polyh_number+1;
1079               int *nodes=getNodesOfPolyhedron(global_polyh_number,nbOfNodes);
1080               for(int i=0;i<nbOfNodes;i++)
1081                 reverse_connectivity[nodes[i]].push_back(global_polyh_number);
1082               delete [] nodes;
1083             }
1084       }
1085
1086     // Full reverse_nodal_connectivity and reverse_nodal_connectivity_index
1087
1088   //calculate size of reverse_nodal_connectivity array
1089     int size_reverse_nodal_connectivity  = 0;
1090     for (int i = 1; i < _numberOfNodes+1; i++)
1091       size_reverse_nodal_connectivity += reverse_connectivity[i].size();
1092
1093     int * reverse_nodal_connectivity_index = new int[_numberOfNodes+1];
1094     int * reverse_nodal_connectivity = new int[size_reverse_nodal_connectivity];
1095
1096     reverse_nodal_connectivity_index[0] = 1;
1097     for (int i = 1; i < _numberOfNodes+1; i++)
1098       {
1099         int size = reverse_connectivity[i].size();
1100         reverse_nodal_connectivity_index[i] = reverse_nodal_connectivity_index[i-1] + size;
1101         for (int j = 0; j < size; j++)
1102           reverse_nodal_connectivity[reverse_nodal_connectivity_index[i-1]-1+j]= reverse_connectivity[i][j];
1103       }
1104
1105     _reverseNodalConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity,
1106                                                     reverse_nodal_connectivity_index,
1107                                                     reverse_nodal_connectivity,true);
1108   }
1109   END_OF(LOC);
1110 }
1111
1112 /*! If not yet done, calculate the Descending Connectivity */
1113 //-------------------------------------------------//
1114 void CONNECTIVITY::calculateDescendingConnectivity()
1115   //-------------------------------------------------//
1116   {
1117   const char * LOC = "CONNECTIVITY::calculateDescendingConnectivity() : ";
1118   BEGIN_OF(LOC);
1119
1120   if (_descending==NULL && _polygonsDescending==NULL && _polyhedronDescending==NULL)
1121     {
1122       if (_nodal==NULL)
1123         {
1124           MESSAGE(LOC<<"No connectivity found !");
1125           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
1126         }
1127       // calcul _descending from _nodal
1128       // we need CONNECTIVITY for constituent
1129
1130       if (_constituent != NULL)
1131         //      throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"ERROR : No descending connectivity defined, but constituent exist !"));
1132         MESSAGE(LOC<<": No descending connectivity defined, but constituent exist !");
1133
1134       if (_entityDimension == 3)
1135         _constituent = new CONNECTIVITY(MED_FACE);
1136       else if (_entityDimension == 2)
1137         _constituent = new CONNECTIVITY(MED_EDGE);
1138       else
1139         {
1140           MESSAGE(LOC<<"We are in 1D");
1141           return;
1142         }
1143
1144       _constituent->_typeConnectivity = MED_NODAL;
1145       _constituent->_numberOfNodes = _numberOfNodes;
1146       // foreach cells, we built array of constituent
1147       int DescendingSize = 0;
1148       for(int i=0; i<_numberOfTypes; i++)
1149         DescendingSize+=(_count[i+1]-_count[i])*_type[i].getNumberOfConstituents(1);
1150       //_descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,DescendingSize);
1151       //const int * descend_connectivity = _descending->getValue();
1152       int * descend_connectivity = new int[DescendingSize];
1153       for (int i=0; i<DescendingSize; i++)
1154         descend_connectivity[i]=0;
1155       //const int * descend_connectivity_index = _descending->getIndex();
1156       int * descend_connectivity_index = new int[_count[_numberOfTypes]];
1157       descend_connectivity_index[0]=1;
1158       medGeometryElement* ConstituentsTypes = new medGeometryElement[2];
1159       ConstituentsTypes[0]=MED_NONE;
1160       ConstituentsTypes[1]=MED_NONE;
1161       int * NumberOfConstituentsForeachType = new int[2];
1162       NumberOfConstituentsForeachType[0]=0;
1163       NumberOfConstituentsForeachType[1]=0;
1164       for(int i=0; i<_numberOfTypes; i++)
1165         {
1166           // initialize descend_connectivity_index array :
1167           int NumberOfConstituents = _type[i].getNumberOfConstituents(1);
1168           for (int j=_count[i];j<_count[i+1];j++)
1169             {
1170               descend_connectivity_index[j]=descend_connectivity_index[j-1]+NumberOfConstituents;
1171               // compute number of constituent of all cell for each type
1172               for(int k=1;k<NumberOfConstituents+1;k++)
1173                 {
1174                   medGeometryElement MEDType = _type[i].getConstituentType(1,k);
1175                   if (ConstituentsTypes[0]==MED_NONE)
1176                     {
1177                       ConstituentsTypes[0]=MEDType;
1178                       NumberOfConstituentsForeachType[0]++;
1179                     } else if (ConstituentsTypes[0]==MEDType)
1180                       NumberOfConstituentsForeachType[0]++;
1181                   else if (ConstituentsTypes[1]==MED_NONE)
1182                     {
1183                       ConstituentsTypes[1]=MEDType;
1184                       NumberOfConstituentsForeachType[1]++;
1185                     } else if (ConstituentsTypes[1]==MEDType)
1186                       NumberOfConstituentsForeachType[1]++;
1187                   else
1188                     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<MEDType<<" is different of "<<ConstituentsTypes[0]<<" and "<<ConstituentsTypes[1]));
1189                 }
1190             }
1191         }
1192
1193       // we could built _constituent !
1194       int TotalNumberOfConstituents = NumberOfConstituentsForeachType[0]+NumberOfConstituentsForeachType[1];
1195       int TotalNumberOfNodes = NumberOfConstituentsForeachType[0]*(ConstituentsTypes[0]%100)+NumberOfConstituentsForeachType[1]*(ConstituentsTypes[1]%100);
1196
1197       //_constituent->_nodal = new MEDSKYLINEARRAY(TotalNumberOfConstituents,TotalNumberOfNodes);
1198
1199       // we use _constituent->_nodal
1200       //const int * ConstituentNodalConnectivity = _constituent->_nodal->getValue();
1201       //const int * ConstituentNodalConnectivityIndex = _constituent->_nodal->getIndex();
1202       int * ConstituentNodalConnectivity = new int[TotalNumberOfNodes];
1203       int * ConstituentNodalConnectivityIndex = new int[TotalNumberOfConstituents+1];
1204       ConstituentNodalConnectivityIndex[0]=1;
1205
1206       _constituent->_entityDimension = _entityDimension-1;
1207       if (ConstituentsTypes[1]==MED_NONE)
1208         _constituent->_numberOfTypes = 1;
1209       else
1210         _constituent->_numberOfTypes = 2;
1211       _constituent->_geometricTypes = new medGeometryElement[_constituent->_numberOfTypes];
1212       _constituent->_type = new CELLMODEL[_constituent->_numberOfTypes];
1213       _constituent->_count = new med_int[_constituent->_numberOfTypes+1];
1214       _constituent->_count[0]=1;
1215       int* tmp_NumberOfConstituentsForeachType = new med_int[_constituent->_numberOfTypes+1];
1216       tmp_NumberOfConstituentsForeachType[0]=0; // to count constituent of each type
1217       for (int i=0; i<_constituent->_numberOfTypes;i++) {
1218         _constituent->_geometricTypes[i]=ConstituentsTypes[i];
1219         _constituent->_count[i+1]=_constituent->_count[i]+NumberOfConstituentsForeachType[i];
1220         CELLMODEL Type(ConstituentsTypes[i]);
1221         _constituent->_type[i]=Type;
1222         tmp_NumberOfConstituentsForeachType[i+1]=NumberOfConstituentsForeachType[i];
1223         for (int j=tmp_NumberOfConstituentsForeachType[i]; j<tmp_NumberOfConstituentsForeachType[i+1]+tmp_NumberOfConstituentsForeachType[i]; j++)
1224           ConstituentNodalConnectivityIndex[j+1]=ConstituentNodalConnectivityIndex[j]+(ConstituentsTypes[i]%100);
1225       }
1226       delete [] ConstituentsTypes;
1227       delete [] NumberOfConstituentsForeachType;
1228
1229       // we need reverse nodal connectivity
1230       if (! _reverseNodalConnectivity)
1231         calculateReverseNodalConnectivity();
1232       const int * ReverseNodalConnectivityIndex = _reverseNodalConnectivity->getIndex();
1233       const int * ReverseNodalConnectivityValue = _reverseNodalConnectivity->getValue();
1234
1235       // array to keep reverse descending connectivity
1236       int * ReverseDescendingConnectivityValue = new int[TotalNumberOfConstituents*2];
1237
1238       int TotalNumberOfSubCell = 0;
1239       for (int i=0; i<_numberOfTypes; i++)
1240         { // loop on all cell type
1241           CELLMODEL Type = _type[i];
1242           //      int NumberOfNodesPerCell = Type.getNumberOfNodes();
1243           int NumberOfConstituentPerCell = Type.getNumberOfConstituents(1);
1244           for (int j=_count[i]; j<_count[i+1]; j++) // we loop on all cell of this type
1245             for (int k=1; k<=NumberOfConstituentPerCell; k++)
1246               { // we loop on all sub cell of it
1247                 if (descend_connectivity[descend_connectivity_index[j-1]+k-2]==0) {
1248                   // it is a new sub cell !
1249                   //        TotalNumberOfSubCell++;
1250                   // Which type ?
1251                   if (Type.getConstituentType(1,k)==_constituent->_geometricTypes[0]){
1252                     tmp_NumberOfConstituentsForeachType[0]++;
1253                     TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[0];
1254                   } else {
1255                     tmp_NumberOfConstituentsForeachType[1]++;
1256                     TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[1];
1257                   }
1258                   //we have maximum two types
1259
1260                   descend_connectivity[descend_connectivity_index[j-1]+k-2]=TotalNumberOfSubCell;
1261                   ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2]=j;
1262                   int NumberOfNodesPerConstituent = Type.getConstituentType(1,k)%100;
1263
1264                   int * NodesLists = new int[NumberOfNodesPerConstituent];
1265                   for (int l=0; l<NumberOfNodesPerConstituent; l++) {
1266                     NodesLists[l]=_nodal->getIJ(j,Type.getNodeConstituent(1,k,l+1));
1267                     ConstituentNodalConnectivity[ConstituentNodalConnectivityIndex[TotalNumberOfSubCell-1]-1+l]=NodesLists[l];
1268                   }
1269                   // we use reverse_nodal_connectivity to find the other element which contain this sub cell
1270
1271                   // all elements which contains first node of sub cell :
1272                   int ReverseNodalConnectivityIndex_0 = ReverseNodalConnectivityIndex[NodesLists[0]-1];
1273                   int ReverseNodalConnectivityIndex_1 = getIndexOfEndClassicElementInReverseNodal(ReverseNodalConnectivityValue,ReverseNodalConnectivityIndex,NodesLists[0]-1);
1274                   //ReverseNodalConnectivityIndex[NodesLists[0]];
1275                   int NumberOfCellsInList = ReverseNodalConnectivityIndex_1-ReverseNodalConnectivityIndex_0;
1276
1277                   if (NumberOfCellsInList > 0)
1278                     { // we could have no element !
1279                       int * CellsList = new int[NumberOfCellsInList];
1280                       for (int l=ReverseNodalConnectivityIndex_0; l<ReverseNodalConnectivityIndex_1; l++)
1281                         CellsList[l-ReverseNodalConnectivityIndex_0]=ReverseNodalConnectivityValue[l-1];
1282
1283                       // foreach node in sub cell, we search elements which are in common
1284                       // at the end, we must have only one !
1285
1286                       for (int l=1; l<NumberOfNodesPerConstituent; l++) {
1287                         int NewNumberOfCellsInList = 0;
1288                         int * NewCellsList = new int[NumberOfCellsInList];
1289                         for (int l1=0; l1<NumberOfCellsInList; l1++)
1290                           for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<getIndexOfEndClassicElementInReverseNodal(ReverseNodalConnectivityValue,ReverseNodalConnectivityIndex,NodesLists[l]-1); l2++)
1291                           //for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<ReverseNodalConnectivityIndex[NodesLists[l]]; l2++)
1292                             {
1293                               if (CellsList[l1]<ReverseNodalConnectivityValue[l2-1])
1294                                 // increasing order : CellsList[l1] are not in elements list of node l
1295                                 break;
1296                               if ((CellsList[l1]==ReverseNodalConnectivityValue[l2-1])&(CellsList[l1]!=j)) {
1297                                 // we have found one
1298                                 NewCellsList[NewNumberOfCellsInList]=CellsList[l1];
1299                                 NewNumberOfCellsInList++;
1300                                 break;
1301                               }
1302                             }
1303                         NumberOfCellsInList = NewNumberOfCellsInList;
1304
1305                         delete [] CellsList;
1306                         CellsList = NewCellsList;
1307                       }
1308
1309                       if (NumberOfCellsInList > 0) { // We have found some elements !
1310                         int CellNumber = CellsList[0];
1311                         //delete [] CellsList;
1312                         if (NumberOfCellsInList>1)
1313                           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one other Cell ("<<NumberOfCellsInList<<") !"));
1314
1315                         if (NumberOfCellsInList==1)
1316                           {
1317                             ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=CellNumber;
1318
1319                             // we search sub cell number in this cell to not calculate it another time
1320                             // which type ?
1321                             CELLMODEL Type2;
1322                             for (int l=0; l<_numberOfTypes; l++)
1323                               if (CellNumber < _count[l+1]) {
1324                                 Type2=_type[l];
1325                                 break;
1326                               }
1327                             //int sub_cell_count2 = Type2.get_entities_count(1);
1328                             //int nodes_cell_count2 = Type2.get_nodes_count();
1329                             bool find2 = false;
1330                             for (int l=1; l<=Type2.getNumberOfConstituents(1);l++) { // on all sub cell
1331                               int counter = 0;
1332                               for (int m=1; m<=Type2.getConstituentType(1,l)%100; m++)
1333                                 for (int n=1; n<=Type.getConstituentType(1,k)%100; n++)
1334                                   {
1335                                     if (_nodal->getIJ(CellNumber,Type2.getNodeConstituent(1,l,m)) == NodesLists[n-1])
1336                                       counter++;
1337                                   }
1338                               if (counter==Type.getConstituentType(1,k)%100)
1339                                 {
1340                                   descend_connectivity[descend_connectivity_index[CellNumber-1]+l-2]=-1*TotalNumberOfSubCell; // because, see it in other side !
1341                                   find2 = true;
1342                                 }
1343                               if (find2)
1344                                 break;
1345                             }
1346                             if (!find2)
1347                               MESSAGE(LOC<<"ERROR ERROR ERROR ERROR ERROR : we find any subcell !!!"); // exception ?
1348                           }
1349                       } else {
1350                         ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=0;
1351                       }
1352                       delete [] CellsList;
1353                     }
1354
1355                   delete [] NodesLists;
1356                 }
1357               }
1358         }
1359       // we adjust _constituent
1360       int NumberOfConstituent=0;
1361       int SizeOfConstituentNodal=0;
1362       for (int i=0;i<_constituent->_numberOfTypes; i++) {
1363         NumberOfConstituent += tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1;
1364         SizeOfConstituentNodal += (tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1)*_constituent->_type[i].getNumberOfNodes();
1365       }
1366       // we built new _nodal attribute in _constituent
1367       //MEDSKYLINEARRAY *ConstituentNodal = new MEDSKYLINEARRAY(NumberOfConstituent,SizeOfConstituentNodal);
1368       //const int *ConstituentNodalValue = ConstituentNodal->getValue();
1369       //const int *ConstituentNodalIndex = ConstituentNodal->getIndex();
1370       int *ConstituentNodalValue = new int[SizeOfConstituentNodal];
1371       int *ConstituentNodalIndex = new int[NumberOfConstituent+1];
1372       ConstituentNodalIndex[0]=1;
1373       // we build _reverseDescendingConnectivity
1374       //_reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,2*NumberOfConstituent);
1375       //const int *reverseDescendingConnectivityValue = _reverseDescendingConnectivity->getValue();
1376       //const int *reverseDescendingConnectivityIndex = _reverseDescendingConnectivity->getIndex();
1377       int *reverseDescendingConnectivityValue = new int[2*NumberOfConstituent];
1378       int *reverseDescendingConnectivityIndex = new int[NumberOfConstituent+1];
1379       reverseDescendingConnectivityIndex[0]=1;
1380
1381       // first constituent type
1382       for(int j=0; j<tmp_NumberOfConstituentsForeachType[0]; j++) {
1383         ConstituentNodalIndex[j+1]=ConstituentNodalConnectivityIndex[j+1];
1384         for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1385           ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[k-1];
1386         }
1387         reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1388         for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1389           reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[k-1];
1390         }
1391       }
1392       // second type if any
1393       if (_constituent->_numberOfTypes==2) {
1394         int offset = _constituent->_count[1]-tmp_NumberOfConstituentsForeachType[0]-1;
1395         int offset1=offset*_constituent->_type[0].getNumberOfNodes();
1396         int offset2=offset*2;
1397         int NumberOfNodesPerConstituent = _constituent->_type[1].getNumberOfNodes();
1398         for(int j=tmp_NumberOfConstituentsForeachType[0]; j<NumberOfConstituent; j++) {
1399           ConstituentNodalIndex[j+1]=ConstituentNodalIndex[j]+NumberOfNodesPerConstituent;
1400           for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1401             ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[offset1+k-1];
1402           }
1403           reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1404           for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1405             reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[offset2+k-1];
1406           }
1407         }
1408         _constituent->_count[2]=NumberOfConstituent+1;
1409         // we correct _descending to adjust face number
1410         for(int j=0;j<DescendingSize;j++)
1411           if (descend_connectivity[j]>tmp_NumberOfConstituentsForeachType[0])
1412             descend_connectivity[j]-=offset;
1413       }
1414
1415       delete [] ConstituentNodalConnectivityIndex;
1416       delete [] ConstituentNodalConnectivity;
1417       delete [] ReverseDescendingConnectivityValue;
1418       _constituent->_count[1]=tmp_NumberOfConstituentsForeachType[0]+1;
1419       delete [] tmp_NumberOfConstituentsForeachType;
1420
1421       _descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,
1422                                         DescendingSize,
1423                                         descend_connectivity_index,
1424                                         descend_connectivity);
1425       delete [] descend_connectivity_index;
1426       delete [] descend_connectivity;
1427
1428       ////
1429       vector<int> PolyDescending;
1430       vector<int> Reversedescendingconnectivityvalue(reverseDescendingConnectivityValue,reverseDescendingConnectivityValue + 2*NumberOfConstituent);
1431       vector<int> Reversedescendingconnectivityindex(reverseDescendingConnectivityIndex,reverseDescendingConnectivityIndex + NumberOfConstituent);
1432       delete [] reverseDescendingConnectivityValue;
1433       delete [] reverseDescendingConnectivityIndex;
1434
1435
1436       // polygons (2D mesh)
1437
1438       vector<int> Constituentnodalvalue(ConstituentNodalValue,ConstituentNodalValue + SizeOfConstituentNodal);
1439       vector<int> Constituentnodalindex(ConstituentNodalIndex,ConstituentNodalIndex + NumberOfConstituent+1);
1440       delete [] ConstituentNodalValue;
1441       delete [] ConstituentNodalIndex;
1442       int NumberOfNewSeg = 0;
1443
1444       for (int i=0; i <getNumberOfPolygons(); i++) // for each polygon
1445         {
1446           const int * vector_begin = &_polygonsNodal->getValue()[_polygonsNodal->getIndex()[i]-1];
1447           int vector_size = _polygonsNodal->getIndex()[i+1]-_polygonsNodal->getIndex()[i]+1;
1448           vector<int> myPolygon(vector_begin, vector_begin+vector_size);
1449           myPolygon[myPolygon.size()-1] = myPolygon[0]; // because first and last point make a segment
1450
1451           for (int j=0; j<myPolygon.size()-1; j++) // for each segment of polygon
1452             {
1453               MEDMODULUSARRAY segment_poly(2,&myPolygon[j]);
1454               int ret_compare = 0;
1455
1456               // we search it in existing segments
1457
1458               for (int k=0; k<Constituentnodalindex.size()-1; k++)
1459                 {
1460                   MEDMODULUSARRAY segment(2,&Constituentnodalvalue[0] + Constituentnodalindex[k]-1);
1461                   ret_compare = segment_poly.compare(segment);
1462                   if (ret_compare) // segment_poly already exists
1463                     {
1464                       PolyDescending.push_back(ret_compare*(k+1)); // we had it to the connectivity
1465                       insert_vector(Reversedescendingconnectivityvalue, 2*k+1, i+1 + getNumberOf(MED_CELL,MED_ALL_ELEMENTS)); // add polygon i to reverse descending connectivity for segment_poly (in 2sd place)
1466                       break;
1467                     }
1468                 }
1469
1470               // segment_poly must be created
1471
1472               if (!ret_compare)
1473                 {
1474                   NumberOfNewSeg++;
1475                   PolyDescending.push_back(NumberOfConstituent+NumberOfNewSeg); // we had it to the connectivity
1476                   Constituentnodalvalue.push_back(segment_poly[0]);
1477                   Constituentnodalvalue.push_back(segment_poly[1]);
1478                   Constituentnodalindex.push_back(Constituentnodalindex[Constituentnodalindex.size()-1] + 2); // we have only segments
1479                   insert_vector(Reversedescendingconnectivityvalue, 2*(NumberOfConstituent+NumberOfNewSeg-1), i+1 + getNumberOf(MED_CELL,MED_ALL_ELEMENTS)); // add polygon i to reverse descending connectivity for segment_poly (in 1st place)
1480                   insert_vector(Reversedescendingconnectivityindex, NumberOfConstituent+NumberOfNewSeg-1, 2*(NumberOfConstituent+NumberOfNewSeg-1)+1); // idem with index
1481                 }
1482             }
1483         }
1484
1485       if (getNumberOfPolygons() > 0)
1486         {
1487           _polygonsDescending = new MEDSKYLINEARRAY(getNumberOfPolygons(),_polygonsNodal->getLength(),_polygonsNodal->getIndex(),&PolyDescending[0]); // index are the same for polygons nodal and descending connectivities
1488
1489           NumberOfConstituent += NumberOfNewSeg;
1490           SizeOfConstituentNodal += 2*NumberOfNewSeg;
1491           _constituent->_count[1] = NumberOfConstituent+1;
1492         }
1493
1494
1495       // polyhedron (3D mesh)
1496
1497       vector<int> Constituentpolygonsnodalvalue;
1498       vector<int> Constituentpolygonsnodalindex(1,1);
1499       int NumberOfNewFaces = 0; // by convention new faces are polygons
1500
1501       for (int i=0; i<getNumberOfPolyhedron(); i++) // for each polyhedron
1502         {
1503           // we create a POLYHEDRONARRAY containing only this polyhedra
1504           int myNumberOfFaces = _polyhedronNodal->getPolyhedronIndex()[i+1]-_polyhedronNodal->getPolyhedronIndex()[i];
1505           int myNumberOfNodes = _polyhedronNodal->getFacesIndex()[_polyhedronNodal->getPolyhedronIndex()[i+1]-1]-_polyhedronNodal->getFacesIndex()[_polyhedronNodal->getPolyhedronIndex()[i]-1];
1506           POLYHEDRONARRAY myPolyhedra(1,myNumberOfFaces,myNumberOfNodes);
1507           vector<int> myFacesIndex(_polyhedronNodal->getFacesIndex() + _polyhedronNodal->getPolyhedronIndex()[i]-1, _polyhedronNodal->getFacesIndex() + _polyhedronNodal->getPolyhedronIndex()[i]-1 + myNumberOfFaces+1);
1508           for (int j=myNumberOfFaces; j>=0; j--)
1509             myFacesIndex[j] -= myFacesIndex[0]-1;
1510           myPolyhedra.setFacesIndex(&myFacesIndex[0]);
1511           myPolyhedra.setNodes(_polyhedronNodal->getNodes() + _polyhedronNodal->getFacesIndex()[_polyhedronNodal->getPolyhedronIndex()[i]-1]-1);
1512
1513           for (int j=0; j<myPolyhedra.getNumberOfFaces(); j++) // for each face of polyhedra
1514             {
1515               int myFaceNumberOfNodes = myPolyhedra.getFacesIndex()[j+1]-myPolyhedra.getFacesIndex()[j];
1516               MEDMODULUSARRAY face_poly(myFaceNumberOfNodes,myPolyhedra.getNodes() + myPolyhedra.getFacesIndex()[j]-1);
1517               int ret_compare = 0;
1518
1519               // we search it in existing faces
1520
1521               // we search first in TRIA3 and QUAD4
1522               if (myFaceNumberOfNodes == 3 || myFaceNumberOfNodes == 4)
1523                 {
1524                   int Begin = -1; // first TRIA3 or QUAD4
1525                   int Number = 0; // number of TRIA3 or QUAD4
1526                   for (int k=0; k<Constituentnodalindex.size()-1; k++)
1527                     if (Constituentnodalindex[k+1]-Constituentnodalindex[k] == myFaceNumberOfNodes)
1528                       {
1529                         if (Begin == -1)
1530                           Begin = k;
1531                         Number++;
1532                       }
1533
1534                   for (int k=0; k<Number; k++)
1535                     {
1536                       MEDMODULUSARRAY face(myFaceNumberOfNodes,&Constituentnodalvalue[0] + Constituentnodalindex[Begin+k]-1);
1537                       ret_compare = face_poly.compare(face);
1538                       if (ret_compare)
1539                         {
1540                           PolyDescending.push_back(ret_compare*(Begin+k+1)); // we had it to the connectivity
1541                           insert_vector(Reversedescendingconnectivityvalue, 2*(Begin+k)+1, i+1 + getNumberOf(MED_CELL,MED_ALL_ELEMENTS)); // add polyhedra i to reverse descending connectivity for face_poly (in 2sd place)
1542                           break;
1543                         }
1544                     }
1545                 }
1546
1547               // we search last in POLYGONS
1548               if (!ret_compare)
1549                 {
1550                   for (int k=0; k<static_cast<int>(Constituentpolygonsnodalindex.size())-1; k++) // we must cast the unsigned int into int before doing -1
1551                     {
1552                       if (Constituentpolygonsnodalindex[k+1]-Constituentpolygonsnodalindex[k] == myFaceNumberOfNodes)
1553                         {
1554                           MEDMODULUSARRAY face(myFaceNumberOfNodes,&Constituentpolygonsnodalvalue[0] + Constituentpolygonsnodalindex[k]-1);
1555                           ret_compare = face_poly.compare(face);
1556                           if (ret_compare)
1557                             {
1558                               PolyDescending.push_back(ret_compare*(NumberOfConstituent+k+1)); // we had it to the connectivity
1559                               insert_vector(Reversedescendingconnectivityvalue, 2*(NumberOfConstituent+k)+1, i+1 + getNumberOf(MED_CELL,MED_ALL_ELEMENTS)); // add polyhedra i to reverse descending connectivity for face_poly (in 2sd place)
1560                               break;
1561                             }
1562                         }
1563                     }
1564                 }
1565
1566               // if not found, face_poly must be created
1567
1568               if (!ret_compare)
1569                 {
1570                   NumberOfNewFaces++;
1571                   PolyDescending.push_back(NumberOfConstituent+NumberOfNewFaces); // we had it to the connectivity
1572                   for (int k=0; k<myFaceNumberOfNodes; k++)
1573                     Constituentpolygonsnodalvalue.push_back(face_poly[k]);
1574                   Constituentpolygonsnodalindex.push_back(Constituentpolygonsnodalindex[Constituentpolygonsnodalindex.size()-1]+myFaceNumberOfNodes);
1575                   insert_vector(Reversedescendingconnectivityvalue, 2*(NumberOfConstituent+NumberOfNewFaces-1), i+1 + getNumberOf(MED_CELL,MED_ALL_ELEMENTS)); // add polyhedra i to reverse descending connectivity for face_poly (in 1st place)
1576                   insert_vector(Reversedescendingconnectivityindex, NumberOfConstituent+NumberOfNewFaces-1, 2*(NumberOfConstituent+NumberOfNewFaces-1)+1); // idem with index
1577                 }
1578             }
1579         }
1580
1581       if (getNumberOfPolyhedron() > 0)
1582         {
1583           _polyhedronDescending = new MEDSKYLINEARRAY(getNumberOfPolyhedron(),_polyhedronNodal->getNumberOfFaces(),_polyhedronNodal->getPolyhedronIndex(),&PolyDescending[0]); // polyhedron index are the same for nodal and descending connectivities
1584
1585           if (_constituent->_polygonsNodal != NULL)
1586             delete [] _constituent->_polygonsNodal;
1587           _constituent->_polygonsNodal = new MEDSKYLINEARRAY(Constituentpolygonsnodalindex.size()-1,Constituentpolygonsnodalvalue.size(),&Constituentpolygonsnodalindex[0],&Constituentpolygonsnodalvalue[0]);
1588         }
1589
1590       // delete _constituent->_nodal
1591       _constituent->_nodal = new MEDSKYLINEARRAY(NumberOfConstituent,
1592                                                  SizeOfConstituentNodal,
1593                                                  &Constituentnodalindex[0],
1594                                                  &Constituentnodalvalue[0]);
1595
1596       int NumberOfConstituentWithPolygons = NumberOfConstituent + NumberOfNewFaces;
1597       Reversedescendingconnectivityindex.push_back(Reversedescendingconnectivityindex[Reversedescendingconnectivityindex.size()-1]+2); // we complete the index
1598       _reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituentWithPolygons+1,
1599                                                            2*NumberOfConstituentWithPolygons,
1600                                                            &Reversedescendingconnectivityindex[0],
1601                                                            &Reversedescendingconnectivityvalue[0]);
1602       ////
1603     }
1604   END_OF(LOC);
1605   }
1606
1607 /*! Not implemented yet */
1608 //--------------------------------------------------------------------//
1609 void CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity)
1610 //--------------------------------------------------------------------//
1611 {
1612   const char * LOC = "CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity) : ";
1613   BEGIN_OF(LOC);
1614
1615   MESSAGE(LOC<<"method not yet implemented " << myConnectivity._entity);
1616   // Mesh dimension !
1617
1618   END_OF(LOC);
1619   return;
1620 }
1621
1622
1623 /*!
1624   Returns the geometry of an element given by its entity type & its global number.
1625
1626   Example : medGeometryElement myType = myConnectivity.getElementType(MED_CELL,35);
1627 */
1628 //--------------------------------------------------------------------//
1629 medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity,int globalNumber) const
1630 //--------------------------------------------------------------------//
1631 {
1632   const char * LOC = "medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity, int globalNumber) const : ";
1633   BEGIN_OF(LOC);
1634   int globalNumberMin = 1;
1635   int globalNumberMax ;
1636
1637   if (_entity==Entity) globalNumberMax = _count[_numberOfTypes];
1638   else if (_constituent!=NULL) globalNumberMax = _constituent->_count[_constituent->_numberOfTypes];
1639   else
1640     throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
1641
1642   // The globalNumber must verify  : 1 <=  globalNumber < _count[_numberOfTypes] (== totalNumberOfElement+1)
1643
1644   if ( (globalNumber < globalNumberMin) || (globalNumber >  globalNumberMax-1 )  )
1645     throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "globalNumber must be between >= |"
1646                                  << globalNumberMin <<"| and <= |" << globalNumberMax-1 ));
1647
1648   if (_entity==Entity) {
1649     for(int i=1; i<=_numberOfTypes;i++)
1650       if (globalNumber<_count[i])
1651         return _geometricTypes[i-1];
1652   }
1653   else if (_constituent!=NULL)
1654     return _constituent->getElementType(Entity,globalNumber);
1655   else
1656     throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
1657   throw MEDEXCEPTION(LOCALIZED("getElementType : Wrong Number !"));
1658
1659   END_OF(LOC);
1660 }
1661
1662 /*!
1663   Idem as getElementType method except that it manages polygon and polyhedron.
1664 */
1665 //--------------------------------------------------------------------//
1666 medGeometryElement CONNECTIVITY::getElementTypeWithPoly(medEntityMesh Entity,int globalNumber) const
1667 //--------------------------------------------------------------------//
1668 {
1669   int globalNumberMaxOfClassicType;
1670   if(_entity==Entity)
1671     {
1672       globalNumberMaxOfClassicType=_count[_numberOfTypes];
1673       if(globalNumber>=1)
1674         {
1675           if(globalNumber<globalNumberMaxOfClassicType)
1676             {
1677               for(int i=1; i<=_numberOfTypes;i++)
1678                 if (globalNumber<_count[i])
1679                   return _geometricTypes[i-1];
1680               throw MEDEXCEPTION("never happens just for compilo");
1681             }
1682           else
1683             {
1684               int localNumberInPolyArray=globalNumber-globalNumberMaxOfClassicType;
1685               int nbOfPol=getNumberOfElementOfPolyType(_entity);
1686               if(localNumberInPolyArray<nbOfPol)
1687                 return getPolyTypeRelativeTo();
1688               else
1689                 throw MEDEXCEPTION("getElementTypeWithPoly : unexisting element type");
1690             }
1691         }
1692       else
1693         throw MEDEXCEPTION("getElementTypeWithPoly : globalNumber < 1");
1694     }
1695   else
1696     {
1697       if(_constituent!=NULL)
1698         return _constituent->getElementTypeWithPoly(Entity,globalNumber);
1699       else
1700         throw MEDEXCEPTION("getElementTypeWithPoly : Entity not defined !");
1701     }
1702 }
1703
1704 ostream & MEDMEM::operator<<(ostream &os, CONNECTIVITY &co)
1705 {
1706     os << endl << "------------- Entity = ";
1707     switch (co._entity)
1708     {
1709         case MED_CELL:
1710             os << "MED_CELL";
1711             break;
1712         case MED_FACE:
1713             os << "MED_FACE";
1714             break;
1715         case MED_EDGE:
1716             os << "MED_EDGE";
1717             break;
1718         case MED_NODE:
1719             os << "MED_NODE";
1720             break;
1721         case MED_ALL_ENTITIES:
1722             os << "MED_ALL_ENTITIES";
1723             break;
1724         default:
1725             os << "Unknown";
1726             break;
1727     }
1728     os  << " -------------\n\nMedConnectivity : ";
1729     switch (co._typeConnectivity)
1730     {
1731         case MED_NODAL:
1732             os << "MED_NODAL\n";
1733             break;
1734         case MED_DESCENDING:
1735             os << "MED_DESCENDING\n";
1736             break;
1737         default:
1738             break;
1739     }
1740     os << "Entity dimension : " << co._entityDimension << endl;
1741     os << "Number of nodes : " << co._numberOfNodes << endl;
1742     os << "Number of types : " << co._numberOfTypes << endl;
1743     for (int i=0; i!=co._numberOfTypes ; ++i)
1744         os << "  -> Type " << co._geometricTypes[i] << " (" << co._type[i].getName() << ") : "
1745            << co._count[i+1]-co._count[i] << " elements" << endl;
1746
1747     if (co._typeConnectivity == MED_NODAL )
1748     {
1749         for (int i=0; i<co._numberOfTypes; i++)
1750         {
1751             os << endl << co._type[i].getName() << " : " << endl;
1752             int numberofelements = co._count[i+1]-co._count[i];
1753             const int * connectivity = co.getConnectivity(co._typeConnectivity, co._entity, co._geometricTypes[i]);
1754             int numberofnodespercell = co._geometricTypes[i]%100;
1755             for (int j=0;j<numberofelements;j++)
1756             {
1757                 os << setw(6) << j+1 << " : ";
1758                 for (int k=0;k<numberofnodespercell;k++)
1759                     os << connectivity[j*numberofnodespercell+k]<<" ";
1760                 os << endl;
1761             }
1762         }
1763
1764         os << endl << "MED_POLYGON : " << co.getNumberOfPolygons() << " polygons" << endl;
1765         if (co.getNumberOfPolygons() > 0)
1766           {
1767             const int* polygonsconnectivity = co.getPolygonsConnectivity(MED_NODAL,co.getEntity());
1768             const int* polygonsconnectivityindex = co.getPolygonsConnectivityIndex(MED_NODAL,co.getEntity());
1769
1770             for (int i=0; i<co.getNumberOfPolygons(); i++)
1771               {
1772                 int numberofnodesperpolygon = polygonsconnectivityindex[i+1]-polygonsconnectivityindex[i];
1773                 for (int j=0; j<numberofnodesperpolygon; j++)
1774                   os << polygonsconnectivity[polygonsconnectivityindex[i]-1+j] << " ";
1775                 os << endl;
1776               }
1777           }
1778
1779         os << endl << "MED_POLYHEDRA : " << co.getNumberOfPolyhedron() << " polyhedron" << endl;
1780         if (co.getNumberOfPolyhedron() > 0)
1781           {
1782             const int* polyhedronconnectivity = co.getPolyhedronConnectivity(MED_NODAL);
1783             const int* polyhedronfacesindex = co.getPolyhedronFacesIndex();
1784             const int* polyhedronindex = co.getPolyhedronIndex(MED_NODAL);
1785
1786             for (int i=0; i<co.getNumberOfPolyhedron(); i++)
1787               {
1788                 int numberoffacesperpolyhedra = polyhedronindex[i+1]-polyhedronindex[i];
1789                 for (int j=0; j<numberoffacesperpolyhedra; j++)
1790                   {
1791                     int numberofnodesperface = polyhedronfacesindex[polyhedronindex[i]-1+j+1]-polyhedronfacesindex[polyhedronindex[i]-1+j];
1792                     for (int k=0; k<numberofnodesperface; k++)
1793                       os << polyhedronconnectivity[polyhedronfacesindex[polyhedronindex[i]-1+j]-1+k] << " ";
1794                     if (j != numberoffacesperpolyhedra-1) os << "| ";
1795                   }
1796                 os << endl;
1797               }
1798           }
1799
1800     }
1801     else if (co._typeConnectivity == MED_DESCENDING)
1802     {
1803         int numberofelements = co.getNumberOf( co._entity , MED_ALL_ELEMENTS);
1804         if (numberofelements > 0)
1805           {
1806             const int *connectivity =  co.getConnectivity( co._typeConnectivity, co._entity, MED_ALL_ELEMENTS);
1807             const int *connectivity_index =  co.getConnectivityIndex( co._typeConnectivity, co._entity );
1808
1809             for ( int j=0; j!=numberofelements; j++ )
1810               {
1811                 os << "Element " << j+1 << " : ";
1812                 for ( int k=connectivity_index[j]; k!=connectivity_index[j+1]; k++ )
1813                   os << connectivity[k-1] << " ";
1814                 os << endl;
1815               }
1816           }
1817
1818         os << endl << "MED_POLYGON : " << co.getNumberOfPolygons() << " polygons" << endl;
1819         if (co.getNumberOfPolygons() > 0)
1820           {
1821             const int* polygonsconnectivity = co.getPolygonsConnectivity(MED_DESCENDING,co.getEntity());
1822             const int* polygonsconnectivityindex = co.getPolygonsConnectivityIndex(MED_DESCENDING,co.getEntity());
1823
1824             for (int i=0; i<co.getNumberOfPolygons(); i++)
1825               {
1826                 int numberofedgesperpolygon = polygonsconnectivityindex[i+1]-polygonsconnectivityindex[i];
1827                 for (int j=0; j<numberofedgesperpolygon; j++)
1828                   os << polygonsconnectivity[polygonsconnectivityindex[i]-1+j] << " ";
1829                 os << endl;
1830               }
1831           }
1832
1833         os << endl << "MED_POLYHEDRA : " << co.getNumberOfPolyhedron() << " polyhedron" << endl;
1834         if (co.getNumberOfPolyhedron() > 0)
1835           {
1836             const int* polyhedronconnectivity = co.getPolyhedronConnectivity(MED_DESCENDING);
1837             const int* polyhedronindex = co.getPolyhedronIndex(MED_DESCENDING);
1838
1839             for (int i=0; i<co.getNumberOfPolyhedron(); i++)
1840               {
1841                 int numberoffacesperpolyhedra = polyhedronindex[i+1]-polyhedronindex[i];
1842                 for (int j=0; j<numberoffacesperpolyhedra; j++)
1843                   os << polyhedronconnectivity[polyhedronindex[i]-1+j] << " ";
1844                 os << endl;
1845               }
1846           }
1847
1848     }
1849
1850     if (co._constituent)
1851         os << endl << *co._constituent << endl;
1852
1853     return os;
1854 }
1855
1856 /*!
1857   method that adds to vector 'nodes' all the nodes of polyhedron with id 'polyhedronId'.
1858   WARNING the returned pointer should be deallocated. Returned nodes and polyhedronId are in form [1,...]
1859  */
1860 int *CONNECTIVITY::getNodesOfPolyhedron(int polyhedronId, int& lgthOfTab) const
1861 {
1862   const int *nodes=getPolyhedronConnectivity(MED_EN::MED_NODAL);
1863   const int *faces=getPolyhedronFacesIndex();
1864   const int *glob=getPolyhedronIndex(MED_EN::MED_NODAL);
1865   int offsetWithClassicType=getNumberOf(_entity,MED_ALL_ELEMENTS);
1866   set<int> retInSet;
1867   int startFace=glob[polyhedronId-offsetWithClassicType-1]-1;
1868   int endFace=glob[polyhedronId-offsetWithClassicType]-1;
1869   int i;
1870   for(i=startFace;i!=endFace;i++)
1871     {
1872       for(int j=faces[i];j<faces[i+1];j++)
1873         retInSet.insert(nodes[j-1]);
1874     }
1875   lgthOfTab=retInSet.size();
1876   int *ret=new int[lgthOfTab];
1877   set<int>::iterator iter=retInSet.begin();
1878   i=0;
1879   for(iter=retInSet.begin();iter!=retInSet.end();iter++)
1880     ret[i++]=*iter;
1881   return ret;
1882 }
1883
1884 /*!
1885   Idem as MESH::getNodesOfPolyhedron except that returned nodes are sorted by face. 'nbOfNodesPerFaces' is an array of size 'nbOfFaces'.
1886   Returned int** has a size of 'nbOfNodesPerFaces' too, and for each element j in int** the size is  nbOfNodesPerFaces[j].
1887   Warning both returned 'nbOfNodesPerFaces' and returned value should be deallocated. Returned nodes and 'polyhedronId' are in form [1,...]
1888  */
1889 int **CONNECTIVITY::getNodesPerFaceOfPolyhedron(int polyhedronId, int& nbOfFaces, int* & nbOfNodesPerFaces) const
1890 {
1891   int i;
1892   const int *nodes=getPolyhedronConnectivity(MED_EN::MED_NODAL);
1893   const int *faces=getPolyhedronFacesIndex();
1894   const int *glob=getPolyhedronIndex(MED_EN::MED_NODAL);
1895   int offsetWithClassicType=getNumberOf(_entity,MED_ALL_ELEMENTS);
1896
1897   int startFace=glob[polyhedronId-offsetWithClassicType-1]-1;
1898   nbOfFaces=glob[polyhedronId-offsetWithClassicType]-startFace-1;
1899   nbOfNodesPerFaces=new int[nbOfFaces];
1900   int **ret=new int* [nbOfFaces];
1901   for(i=0;i<nbOfFaces;i++)
1902     {
1903       int startNode=faces[startFace+i]-1;
1904       nbOfNodesPerFaces[i]=faces[startFace+i+1]-startNode-1;
1905       ret[i]=(int *)(nodes)+startNode;
1906     }
1907   return ret;
1908 }
1909
1910 int CONNECTIVITY::getNumberOfTypesWithPoly(MED_EN::medEntityMesh Entity) const
1911 {
1912   if (_entity==Entity)
1913     return _numberOfTypes+getNumberOfPolyType();
1914   else if (_constituent!=NULL)
1915     return _constituent->getNumberOfTypesWithPoly(Entity);
1916   else
1917     return 0;
1918 }
1919
1920 int CONNECTIVITY::getNumberOfPolyType()  const
1921 {
1922   if (_entity==MED_CELL && _entityDimension==3)
1923     {
1924       if(getNumberOfPolyhedron()>0)
1925         return 1;
1926     }
1927   else if ((_entity==MED_CELL && _entityDimension==2) || (_entity==MED_FACE && _entityDimension==2))
1928     if(getNumberOfPolygons()>0)
1929       return 1;
1930   return 0;
1931 }
1932
1933 int CONNECTIVITY::getNumberOfElementOfPolyType(MED_EN::medEntityMesh Entity)  const
1934 {
1935   if(_entity==Entity)
1936     {
1937       if (_entity==MED_CELL && _entityDimension==3)
1938         return getNumberOfPolyhedron();
1939       else if ((_entity==MED_CELL && _entityDimension==2) || (_entity==MED_FACE && _entityDimension==2))
1940         return getNumberOfPolygons();
1941       return 0;
1942     }
1943   else
1944     {
1945       if(_constituent!=NULL)
1946         return _constituent->getNumberOfElementOfPolyType(Entity);
1947       else
1948         throw MEDEXCEPTION("_constituent required to evaluate getNumberOfElementOfPolyType");
1949     }
1950 }
1951
1952 /*
1953   Method equivalent to getGeometricTypes except that it includes not only classical Types but polygons/polyhedra also.
1954   WARNING the returned array MUST be deallocated.
1955  */
1956 MED_EN::medGeometryElement * CONNECTIVITY::getGeometricTypesWithPoly(MED_EN::medEntityMesh Entity) const throw (MEDEXCEPTION)
1957 {
1958   if(_entity==Entity)
1959     {
1960       int nbOfTypesTotal=getNumberOfTypesWithPoly(Entity);
1961       int nbOfTypesWithoutPoly=getNumberOfTypes(Entity);
1962       medGeometryElement *ret=new medGeometryElement[nbOfTypesTotal];
1963       memcpy(ret,getGeometricTypes(Entity),nbOfTypesWithoutPoly*sizeof(medGeometryElement));
1964       if(nbOfTypesTotal>nbOfTypesWithoutPoly)
1965         {
1966           if (Entity==MED_CELL && _entityDimension==3)
1967             ret[nbOfTypesWithoutPoly]=MED_EN::MED_POLYHEDRA;
1968           else if((Entity==MED_CELL && _entityDimension==2) || (Entity==MED_FACE && _entityDimension==2))
1969             ret[nbOfTypesWithoutPoly]=MED_EN::MED_POLYGON;
1970         }
1971       return ret;
1972     }
1973   else
1974     {
1975       if(_constituent)
1976         return _constituent->getGeometricTypesWithPoly(Entity);
1977       throw MEDEXCEPTION("_constituent required to evaluate getGeometricTypesWithPoly");
1978     }
1979 }
1980
1981 /*
1982   Method used in CalculateDescendingConnectivity. So it's typically a private method.
1983   The aim of this method is to hide to CalculateDescendingConnectivity algorithm the fact that in reverse connectivity polygons and polyhedrons
1984   are not in separate data structure, contrary to not reverse connectivities.
1985  */
1986 int CONNECTIVITY::getIndexOfEndClassicElementInReverseNodal(const int *reverseNodalValue, const int *reverseNodalIndex, int rk)  const
1987 {
1988   int nbOfLastElt=getNumberOf(_entity,MED_ALL_ELEMENTS);
1989   int ret=reverseNodalIndex[rk];
1990   for(int i=reverseNodalIndex[rk];i<reverseNodalIndex[rk+1];i++)
1991     {
1992       if(reverseNodalValue[i-1]<=nbOfLastElt)
1993         ret++;
1994     }
1995   return ret;
1996 }
1997
1998 /*
1999   Method that inverts the face with faceId 'faceId' in the data structure. If it's a polygon face 'faceId' is a value between 1 and nbOfPolygons.
2000   This method has to be applied on a instance of MEDMEM::CONNECTIVITY with _entityDimension==3.
2001   WARNING calculateDescendingConnectivity must have been called before.
2002  */
2003 void CONNECTIVITY::invertConnectivityForAFace(int faceId, const int *nodalConnForFace, bool polygonFace)
2004 {
2005   // we have always 2 neighbourings
2006   int cell1 = _reverseDescendingConnectivity->getIJ(faceId,1);
2007   int cell2 = _reverseDescendingConnectivity->getIJ(faceId,2);
2008
2009   if (cell2 != 0) { // we are not on border, update compulsary. If we aren't on border no update necessary so WARNING because user specified a bad oriented face
2010     // Updating _reverseDescendingConnectivity
2011     _reverseDescendingConnectivity->setIJ(faceId,1,cell2);
2012     _reverseDescendingConnectivity->setIJ(faceId,2,cell1);
2013     // Updating _constituent->_nodal because of reversity
2014     MEDSKYLINEARRAY *currentNodal=(!polygonFace)?_constituent->_nodal:_constituent->_polygonsNodal;
2015     MEDSKYLINEARRAY *currentDescending=(!polygonFace)?_descending:_polygonsDescending;
2016     const int *descendingNodalIndex=(!polygonFace)?_constituent->_nodal->getIndex():_constituent->_polygonsNodal->getIndex();
2017     const int *newDescendingIndex=(!polygonFace)?_descending->getIndex():_polygonsDescending->getIndex();
2018     for(int iarray=1;iarray<=(descendingNodalIndex[faceId]-descendingNodalIndex[faceId-1]);iarray++)
2019       currentNodal->setIJ(faceId,iarray,nodalConnForFace[iarray-1]);
2020
2021     // Updating _descending for cell1 and cell2
2022     for(int iface=newDescendingIndex[cell1-1];iface<=newDescendingIndex[cell1];iface++)
2023       if (currentDescending->getIndexValue(iface)==faceId)
2024         currentDescending->setIndexValue(iface,-faceId);
2025       else if (currentDescending->getIndexValue(iface)==-faceId)
2026         currentDescending->setIndexValue(iface,faceId);
2027
2028     for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];iface++)
2029       if (currentDescending->getIndexValue(iface)==faceId)
2030         currentDescending->setIndexValue(iface,-faceId);
2031       else if (_descending->getIndexValue(iface)==-faceId)
2032         currentDescending->setIndexValue(iface,faceId);
2033   }
2034 }
2035
2036 /*
2037   Method with 2 output : the connectivity required and its length 'lgth'. This method gives the connectivity independently it is a polygons/polyhedrons or normal
2038   element.
2039  */
2040 const int * CONNECTIVITY::getConnectivityOfAnElementWithPoly(MED_EN::medConnectivity ConnectivityType, MED_EN::medEntityMesh Entity, int Number, int &lgth)
2041 {
2042   if(Entity==MED_EN::MED_NODE)
2043     throw  MEDEXCEPTION("No connectivity attached to a node entity");
2044   if(Entity==_entity)
2045     {
2046       if(_entity==MED_EDGE && _entityDimension==1)
2047         {
2048           const int * newConstituentValue = 0;
2049           const int * newConstituentIndex = 0;
2050           newConstituentValue = getConnectivity(ConnectivityType,Entity,MED_ALL_ELEMENTS);
2051           newConstituentIndex = getConnectivityIndex(ConnectivityType,Entity);
2052           lgth=newConstituentIndex[Number]-newConstituentIndex[Number-1];
2053           return newConstituentValue+newConstituentIndex[Number-1]-1;
2054         }
2055       int nbOfClassicalElements=getNumberOf(_entity,MED_ALL_ELEMENTS);
2056       if((_entity==MED_FACE && _entityDimension==2) || (_entity==MED_CELL && _entityDimension==2))
2057         {
2058           const int * newConstituentValue = 0;
2059           const int * newConstituentIndex = 0;
2060           if(Number<=nbOfClassicalElements)
2061             {
2062               newConstituentValue = getConnectivity(ConnectivityType,Entity,MED_ALL_ELEMENTS);
2063               newConstituentIndex = getConnectivityIndex(ConnectivityType,Entity);
2064               lgth=newConstituentIndex[Number]-newConstituentIndex[Number-1];
2065               return newConstituentValue+newConstituentIndex[Number-1]-1;
2066             }
2067           else
2068             {
2069               int localNumber=Number-nbOfClassicalElements-1;
2070               if(localNumber<getNumberOfPolygons())
2071                 {
2072                   newConstituentValue = getPolygonsConnectivity(ConnectivityType,Entity);
2073                   newConstituentIndex = getPolygonsConnectivityIndex(ConnectivityType,Entity);
2074                   lgth=newConstituentIndex[localNumber+1]-newConstituentIndex[localNumber];
2075                   return newConstituentValue+newConstituentIndex[localNumber]-1;
2076                 }
2077               else
2078                 throw  MEDEXCEPTION("Unknown number");
2079             }
2080         }
2081       else if(_entity==MED_CELL && _entityDimension==3)
2082         {
2083           const int * newConstituentValue = 0;
2084           const int * newConstituentIndex = 0;
2085           if(Number<=nbOfClassicalElements)
2086             {
2087               newConstituentValue = getConnectivity(ConnectivityType,Entity,MED_ALL_ELEMENTS);
2088               newConstituentIndex = getConnectivityIndex(ConnectivityType,Entity);
2089               lgth=newConstituentIndex[Number]-newConstituentIndex[Number-1];
2090               return newConstituentValue+newConstituentIndex[Number-1]-1;
2091             }
2092           else
2093             {
2094               int localNumber=Number-nbOfClassicalElements-1;
2095               if(localNumber<getNumberOfPolyhedron())
2096                 {
2097                   if(ConnectivityType==MED_NODAL)
2098                     throw  MEDEXCEPTION("NODAL Connectivity required for a polyhedron");
2099 //                newConstituentValue = _polyhedronDescending->getValue();
2100 //                newConstituentIndex = _polyhedronDescending->getIndex();
2101                   newConstituentValue = getPolyhedronConnectivity( ConnectivityType );
2102                   newConstituentIndex = getPolyhedronIndex( ConnectivityType );
2103                   lgth=newConstituentIndex[localNumber+1]-newConstituentIndex[localNumber];
2104                   return newConstituentValue+newConstituentIndex[localNumber]-1;
2105                 }
2106               else
2107                 throw  MEDEXCEPTION("Unknown number");
2108             }
2109         }
2110       else
2111         throw MEDEXCEPTION("No connectivity available");
2112     }
2113   else
2114     {
2115       if(_constituent==NULL)
2116         calculateDescendingConnectivity();
2117       return _constituent->getConnectivityOfAnElementWithPoly(ConnectivityType,Entity,Number,lgth);
2118     }
2119 }
2120
2121 int CONNECTIVITY::getNumberOfElementsWithPoly(MED_EN::medEntityMesh Entity, MED_EN::medGeometryElement Type) const
2122 {
2123   if(Entity==MED_EN::MED_NODE)
2124     return _numberOfNodes;
2125   if(Entity==_entity)
2126     {
2127       if(Type==MED_POLYGON || Type==MED_POLYHEDRA)
2128         return getNumberOfElementOfPolyType(_entity);
2129       else
2130         return getNumberOf(_entity,Type);
2131     }
2132   else
2133     {
2134       if(_constituent)
2135         return _constituent->getNumberOfElementsWithPoly(Entity,Type);
2136       else
2137         throw MEDEXCEPTION("CONNECTIVITY::getNumberOfElementsWithPoly : _constituent needed");
2138     }
2139 }
2140
2141 /*!
2142   Perform a deep comparison of the 2 connectivities in NODAL mode and on all elements.
2143 */
2144 bool CONNECTIVITY::deepCompare(const CONNECTIVITY& other) const
2145 {
2146   CONNECTIVITY* temp=(CONNECTIVITY* )this;
2147   const int *conn1=temp->getConnectivity(MED_NODAL,_entity,MED_ALL_ELEMENTS);
2148   int size1=temp->getConnectivityLength(MED_NODAL,_entity,MED_ALL_ELEMENTS);
2149   temp=(CONNECTIVITY* )(&other);
2150   const int *conn2=temp->getConnectivity(MED_NODAL,_entity,MED_ALL_ELEMENTS);
2151   int size2=temp->getConnectivityLength(MED_NODAL,_entity,MED_ALL_ELEMENTS);
2152   if(size1!=size2)
2153     return false;
2154   bool ret=true;
2155   for(int i=0;i<size1 && ret;i++)
2156     {
2157       ret=(conn1[i]==conn2[i]);
2158     }
2159   return ret;
2160 }