src/MEDMEM/MEDMEM_Connectivity.cxx

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