src/MEDMEM/MEDMEM_Connectivity.cxx

   1 //  Copyright (C) 2007-2008  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 #include "MEDMEM_Connectivity.hxx"
  23 #include "MEDMEM_Family.hxx"
  24 #include "MEDMEM_Group.hxx"
  25 #include "MEDMEM_CellModel.hxx"
  26
  27 #include "MEDMEM_SkyLineArray.hxx"
  28 #include "MEDMEM_ModulusArray.hxx"
  29
  30 #include "MEDMEM_STRING.hxx"
  31 #include "MEDMEM_Utilities.hxx"
  32 #include <iomanip>
  33
  34 using namespace std;
  35 using namespace MEDMEM;
  36 using namespace MED_EN;
  37
  38 // Enlarge the vector if necessary to insert the element
  39 static inline void insert_vector(vector<int> &Vect, int Indice, int Element)
  40 {
  41   if (Indice >=(int) Vect.capacity())
  42     Vect.reserve(Indice + 1000);
  43
  44   if (Indice >=(int) Vect.size())
  45     Vect.resize(Indice+1);
  46
  47   Vect[Indice] = Element;
  48 }
  49
  50 void mergeOrderedTabs(const int *tab1, int lgth1, const int *tab2, int lgth2, int *result, int& lgth)
  51 {
  52   int cpt[2]={0,0};
  53   lgth=0;
  54   unsigned char switcher=0;
  55   const int *tabS[2]={tab1,tab2};
  56   while(cpt[0]<lgth1 && cpt[1]<lgth2)
  57     {
  58       if(tabS[1-switcher][cpt[1-switcher]]<tabS[switcher][cpt[switcher]])
  59         cpt[1-switcher]++;
  60       else if(tabS[1-switcher][cpt[1-switcher]]>tabS[switcher][cpt[switcher]])
  61         switcher=1-switcher;
  62       else
  63         {
  64           int tmp=tabS[switcher][cpt[switcher]];
  65           cpt[switcher]++; cpt[1-switcher]++;
  66           result[lgth++]=tmp;
  67         }
  68     }
  69 }
  70
  71 /*!
  72    Default Constructor. /n
  73    Default for Entity is MED_CELL and type of Connectivity is MED_NODAL */
  74 //--------------------------------------------------------------//
  75 CONNECTIVITY::CONNECTIVITY(medEntityMesh Entity /* =MED_CELL */) :
  76 //--------------------------------------------------------------//
  77                                 _entity(Entity),
  78                                 _typeConnectivity(MED_NODAL),
  79                                 _numberOfTypes(0),
  80                                 _geometricTypes((medGeometryElement*)NULL),
  81                                 _type((CELLMODEL*)NULL),
  82                                 _entityDimension(0),
  83                                 _numberOfNodes(0),
  84                                 _count((int*)NULL),
  85                                 _nodal((MEDSKYLINEARRAY*)NULL),
  86                                 _polygonsNodal((MEDSKYLINEARRAY*)NULL),
  87                                 _polyhedronNodal((POLYHEDRONARRAY*)NULL),
  88                                 _descending((MEDSKYLINEARRAY*)NULL),
  89                                 _polygonsDescending((MEDSKYLINEARRAY*)NULL),
  90                                 _polyhedronDescending((MEDSKYLINEARRAY*)NULL),
  91                                 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
  92                                 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
  93                                 _neighbourhood((MEDSKYLINEARRAY*)NULL),
  94                                 _constituent((CONNECTIVITY*)NULL),
  95                                 _isDescendingConnectivityPartial(false)
  96 {
  97   const char* LOC = "CONNECTIVITY(medEntityMesh Entity=MED_CELL)";
  98   BEGIN_OF_MED(LOC);
  99    MESSAGE_MED("CONNECTIVITY(medEntityMesh Entity=MED_CELL)");
 100   _count = new int[1];
 101   _count[0]=1;
 102   END_OF_MED(LOC);
 103 }
 104
 105 /*!
 106    Constructor. /n
 107    Default for Entity is MED_CELL */
 108 //------------------------------------------------------------------------------//
 109 CONNECTIVITY::CONNECTIVITY(int numberOfTypes,medEntityMesh Entity /* =MED_CELL */):
 110 //------------------------------------------------------------------------------//
 111                                 _entity(Entity),
 112                                 _typeConnectivity(MED_NODAL),
 113                                 _numberOfTypes(numberOfTypes),
 114                                 _entityDimension(0),
 115                                 _numberOfNodes(0),
 116                                 _nodal((MEDSKYLINEARRAY*)NULL),
 117                                 _polygonsNodal((MEDSKYLINEARRAY*)NULL),
 118                                 _polyhedronNodal((POLYHEDRONARRAY*)NULL),
 119                                 _descending((MEDSKYLINEARRAY*)NULL),
 120                                 _polygonsDescending((MEDSKYLINEARRAY*)NULL),
 121                                 _polyhedronDescending((MEDSKYLINEARRAY*)NULL),
 122                                 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
 123                                 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
 124                                 _neighbourhood((MEDSKYLINEARRAY*)NULL),
 125                                 _constituent((CONNECTIVITY*)NULL),
 126                                 _isDescendingConnectivityPartial(false)
 127 {
 128   MESSAGE_MED("CONNECTIVITY(int numberOfTypes,medEntityMesh Entity=MED_CELL)");
 129   _geometricTypes = new medGeometryElement[numberOfTypes];
 130   _type = new CELLMODEL[numberOfTypes];
 131   _count = new int[numberOfTypes+1];
 132   _count[0]=1;
 133 }
 134
 135 /*!
 136   Copy Constructor.
 137 */
 138 //------------------------------------------------------//
 139 CONNECTIVITY::CONNECTIVITY (const CONNECTIVITY & m):
 140 //----------------------------------------------------//
 141                                 _entity           (m._entity),
 142                                 _typeConnectivity (m._typeConnectivity),
 143                                 _numberOfTypes    (m._numberOfTypes),
 144                                 _entityDimension  (m._entityDimension),
 145                                 _numberOfNodes    (m._numberOfNodes),
 146                                 _isDescendingConnectivityPartial(m._isDescendingConnectivityPartial)
 147 {
 148  if (m._geometricTypes != NULL)
 149  {
 150     _geometricTypes = new medGeometryElement[_numberOfTypes];
 151     memcpy(_geometricTypes,m._geometricTypes,_numberOfTypes*sizeof(medGeometryElement));
 152  }
 153  else
 154     _geometricTypes = (medGeometryElement *) NULL;
 155
 156  if (m._type != NULL)
 157  {
 158     _type = new CELLMODEL[_numberOfTypes];
 159     for (int i=0;i<_numberOfTypes;i++)
 160         _type[i] = CELLMODEL(m._type[i]);
 161  }
 162  else
 163     _type = (CELLMODEL *) NULL;
 164
 165  if (m._count != NULL)
 166  {
 167       _count = new int[_numberOfTypes+1];
 168       memcpy(_count,m._count,(_numberOfTypes+1)*sizeof(int));
 169  }
 170  else
 171     _count = (int *) NULL;
 172
 173  if (m._nodal != NULL)
 174     _nodal = new MEDSKYLINEARRAY(* m._nodal);
 175  else
 176     _nodal = (MEDSKYLINEARRAY *) NULL;
 177
 178  if (m._polygonsNodal != NULL)
 179    _polygonsNodal = new MEDSKYLINEARRAY(* m._polygonsNodal);
 180  else
 181    _polygonsNodal = (MEDSKYLINEARRAY *) NULL;
 182
 183  if (m._polyhedronNodal != NULL)
 184    _polyhedronNodal = new POLYHEDRONARRAY(* m._polyhedronNodal);
 185  else
 186    _polyhedronNodal = (POLYHEDRONARRAY *) NULL;
 187
 188  if (m._descending != NULL)
 189     _descending = new MEDSKYLINEARRAY(* m._descending);
 190  else
 191     _descending = (MEDSKYLINEARRAY *) NULL;
 192
 193  if (m._polygonsDescending != NULL)
 194    _polygonsDescending = new MEDSKYLINEARRAY(* m._polygonsDescending);
 195  else
 196    _polygonsDescending = (MEDSKYLINEARRAY *) NULL;
 197
 198  if (m._polyhedronDescending != NULL)
 199    _polyhedronDescending = new MEDSKYLINEARRAY(* m._polyhedronDescending);
 200  else
 201    _polyhedronDescending = (MEDSKYLINEARRAY *) NULL;
 202
 203  if (m._reverseNodalConnectivity != NULL)
 204     _reverseNodalConnectivity = new MEDSKYLINEARRAY(* m._reverseNodalConnectivity);
 205  else
 206     _reverseNodalConnectivity = (MEDSKYLINEARRAY *) NULL;
 207
 208  if (m._reverseDescendingConnectivity != NULL)
 209     _reverseDescendingConnectivity = new MEDSKYLINEARRAY(* m._reverseDescendingConnectivity);
 210  else
 211     _reverseDescendingConnectivity = (MEDSKYLINEARRAY *) NULL;
 212
 213  if (m._neighbourhood != NULL)
 214     _neighbourhood = new MEDSKYLINEARRAY(* m._neighbourhood);
 215  else
 216     _neighbourhood = (MEDSKYLINEARRAY *) NULL;
 217
 218  if (m._constituent != NULL)
 219     _constituent = new CONNECTIVITY(* m._constituent);
 220  else
 221     _constituent = (CONNECTIVITY *) NULL;
 222 }
 223
 224 /*!
 225    Destructor./n
 226    desallocates existing pointers */
 227 //----------------------------//
 228 CONNECTIVITY::~CONNECTIVITY()
 229 //----------------------------//
 230 {
 231   MESSAGE_MED("Destructeur de CONNECTIVITY()");
 232
 233   if (_geometricTypes != NULL)
 234      delete [] _geometricTypes;
 235   if (_type != NULL)
 236      delete [] _type;
 237   if (_count != NULL)
 238      delete [] _count;
 239   if (_nodal != NULL)
 240      delete _nodal;
 241   if (_polygonsNodal != NULL)
 242      delete _polygonsNodal;
 243   if (_polyhedronNodal != NULL)
 244      delete _polyhedronNodal;
 245   if (_descending != NULL)
 246      delete _descending;
 247   if (_polygonsDescending != NULL)
 248      delete _polygonsDescending;
 249   if (_polyhedronDescending != NULL)
 250     delete _polyhedronDescending;
 251   if (_reverseNodalConnectivity != NULL)
 252      delete _reverseNodalConnectivity;
 253   if (_reverseDescendingConnectivity != NULL)
 254      delete _reverseDescendingConnectivity;
 255   if (_constituent != NULL)
 256      delete _constituent;
 257 }
 258
 259 /*!
 260    set _constituent to Constituent
 261    be aware desallocation of _constituent is done by CONNECTIVITY:~CONNECTIVITY
 262    throws an exception if Constituent = MED_CELL
 263    A DOCUMENTER
 264     */
 265 //----------------------------------------------------------//
 266 void CONNECTIVITY::setConstituent(CONNECTIVITY * Constituent)
 267                     throw (MEDEXCEPTION)
 268 //----------------------------------------------------------//
 269 {
 270   medEntityMesh Entity = Constituent->getEntity();
 271   if (Entity == MED_CELL)
 272     throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : could not set constituent on MED_CELL !"));
 273
 274   if ((Entity == MED_EDGE)&(_entityDimension == 3))
 275   {
 276     if (_constituent == NULL)
 277       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : Entity not found !"));
 278     _constituent->setConstituent(Constituent);
 279   }
 280   else
 281     {
 282       delete _constituent;
 283       _constituent = Constituent;
 284     }
 285 }
 286
 287 /*! Duplicated Types array in CONNECTIVITY object. */
 288 //--------------------------------------------------------------------//
 289 void CONNECTIVITY::setGeometricTypes(const medGeometryElement * Types,
 290                                      const medEntityMesh Entity)
 291                                      throw (MEDEXCEPTION)
 292 //--------------------------------------------------------------------//
 293 {
 294   if (Entity == _entity)
 295     for (int i=0; i<_numberOfTypes; i++)
 296     {
 297       _geometricTypes[i] = Types[i];
 298       _type[i] = CELLMODEL(_geometricTypes[i]);
 299     }
 300   else
 301     {
 302     if (_constituent == NULL)
 303         throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setGeometricTypes : Entity not found !"));
 304     _constituent->setGeometricTypes(Types,Entity);
 305   }
 306 }
 307
 308 /*! A DOCUMENTER */
 309 //--------------------------------------------------------------------//
 310 void CONNECTIVITY::setCount(const int * Count, const medEntityMesh Entity)
 311                             throw (MEDEXCEPTION)
 312 //--------------------------------------------------------------------//
 313 {
 314   if (Entity == _entity)
 315   {
 316     // this exception was added at 1.18.2.2.2.4 for NPAL17043: "Correction of
 317     // MEDMEM CPPUNIT tests. Integrated a work of V. Bergeaud and A. Geay."
 318     // and removed for PAL19744: "regression in MEDMEM_Connectivity.cxx"
 319     // Commenting this exception at least looks safe as the case
 320     // _numberOfTypes==0 is previewed here
 321     //if (_numberOfTypes==0)
 322     //  throw MEDEXCEPTION("Number of Types was not set before setting counts");
 323     int * index = new int[Count[_numberOfTypes]];
 324     index[0]=1;
 325     _count[0]=1;
 326     for (int i=0; i<_numberOfTypes; i++) {
 327       _count[i+1] = Count[i+1];
 328       int NumberOfNodesPerElement = _type[i].getNumberOfNodes() ;
 329       for (int j=_count[i]; j<_count[i+1]; j++)
 330         index[j] = index[j-1]+NumberOfNodesPerElement;
 331     }
 332     // allocate _nodal
 333     if (_nodal != NULL) delete _nodal;
 334     if (_numberOfTypes != 0)
 335     {
 336       _nodal = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,index[_count[_numberOfTypes]-1]-1);
 337       _nodal->setIndex(index);
 338     }
 339     else
 340       _nodal=0;
 341     delete [] index;
 342   }
 343   else
 344   {
 345     if (_constituent == NULL)
 346       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setCount : Entity not found !"));
 347     _constituent->setCount(Count,Entity);
 348   }
 349 }
 350
 351 //--------------------------------------------------------------------//
 352 void CONNECTIVITY::setNodal(const int * Connectivity,
 353                             const medEntityMesh Entity,
 354                             const medGeometryElement Type)
 355                             throw (MEDEXCEPTION)
 356 //--------------------------------------------------------------------//
 357 {
 358   if (_entity == Entity)
 359   {
 360     // find geometric type
 361     bool find = false;
 362     for (int i=0; i<_numberOfTypes; i++)
 363       if (_geometricTypes[i] == Type) {
 364         find = true;
 365         int NumberOfNodePerElements = _type[i].getNumberOfNodes() ;
 366         //_nodal->setI(i+1,Connectivity);
 367         for( int j=_count[i];j<_count[i+1]; j++)
 368           _nodal->setI(j,Connectivity+(j-_count[i])*NumberOfNodePerElements);
 369       }
 370     if (!find)
 371       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : geometric type not found"));
 372   } else
 373   {
 374     if (_constituent == NULL)
 375       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : Entity not found !"));
 376       _constituent->setNodal(Connectivity,Entity,Type);
 377   }
 378 }
 379
 380
 381 //--------------------------------------------------------------------//
 382 void CONNECTIVITY::setPolygonsConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, const int* PolygonsConnectivity, const int* PolygonsConnectivityIndex, int ConnectivitySize, int NumberOfPolygons)
 383 //--------------------------------------------------------------------//
 384 {
 385   const char* LOC = "CONNECTIVITY::setPolygonsConnectivity";
 386   BEGIN_OF_MED(LOC);
 387
 388   if (_entity == Entity)
 389     {
 390       MEDSKYLINEARRAY* Connectivity = new MEDSKYLINEARRAY(NumberOfPolygons,ConnectivitySize,PolygonsConnectivityIndex,PolygonsConnectivity);
 391
 392       if (ConnectivityType == MED_NODAL)
 393         {
 394           if (_polygonsNodal != (MEDSKYLINEARRAY*) NULL)
 395             delete _polygonsNodal;
 396           _polygonsNodal = Connectivity;
 397         }
 398       else
 399         {
 400           if (_typeConnectivity != MED_DESCENDING)
 401             _typeConnectivity = MED_DESCENDING; //by default it is set to MED_NODAL
 402           if (_polygonsDescending != (MEDSKYLINEARRAY*) NULL)
 403             delete _polygonsDescending;
 404           _polygonsDescending = Connectivity;
 405         }
 406     }
 407   else
 408     {
 409       if (_entityDimension==3)
 410       {
 411         if (_constituent == (CONNECTIVITY*) NULL)
 412           _constituent=new CONNECTIVITY(MED_EN::MED_FACE);
 413       }
 414       else if (_constituent == (CONNECTIVITY*) NULL)
 415       {
 416         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not found !"));
 417       }
 418       _constituent->setPolygonsConnectivity(ConnectivityType, Entity,
 419                                             PolygonsConnectivity, PolygonsConnectivityIndex,
 420                                             ConnectivitySize, NumberOfPolygons);
 421     }
 422 }
 423
 424
 425 //--------------------------------------------------------------------//
 426 void CONNECTIVITY::setPolyhedronConnectivity(medConnectivity ConnectivityType, const int* PolyhedronConnectivity, const int* PolyhedronIndex, int ConnectivitySize, int NumberOfPolyhedron, const int* PolyhedronFacesIndex /* =(int*)NULL */, int NumberOfFaces /* =0 */)
 427 //--------------------------------------------------------------------//
 428 {
 429   const char* LOC = "CONNECTIVITY::setPolyhedronConnectivity";
 430   BEGIN_OF_MED(LOC);
 431
 432   if (_entity == MED_CELL)
 433     {
 434       if (ConnectivityType == MED_NODAL)
 435         {
 436           if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
 437             delete _polyhedronNodal;
 438           _polyhedronNodal = new POLYHEDRONARRAY(NumberOfPolyhedron,NumberOfFaces,ConnectivitySize);
 439 #if defined(IRIX64) || defined(OSF1) || defined(VPP5000) || defined(PCLINUX64)
 440           int i ;
 441           MED_EN::med_int * tmp_PolyhedronIndex = new med_int[NumberOfPolyhedron+1] ;
 442           for ( i = 0 ; i < NumberOfPolyhedron+1 ; i++ )
 443              tmp_PolyhedronIndex[i] = PolyhedronIndex[i] ;
 444           _polyhedronNodal->setPolyhedronIndex(tmp_PolyhedronIndex);
 445           delete [] tmp_PolyhedronIndex ;
 446           MED_EN::med_int * tmp_PolyhedronFacesIndex = new med_int[NumberOfFaces+1] ;
 447           for ( i = 0 ; i < NumberOfFaces+1 ; i++ )
 448              tmp_PolyhedronFacesIndex[i] = PolyhedronFacesIndex[i] ;
 449           _polyhedronNodal->setFacesIndex(tmp_PolyhedronFacesIndex);
 450           delete [] tmp_PolyhedronFacesIndex ;
 451           MED_EN::med_int * tmp_PolyhedronConnectivity = new med_int[ConnectivitySize] ;
 452           for ( i = 0 ; i < ConnectivitySize ; i++ )
 453              tmp_PolyhedronConnectivity[i] = PolyhedronConnectivity[i] ;
 454           _polyhedronNodal->setNodes(tmp_PolyhedronConnectivity);
 455           delete [] tmp_PolyhedronConnectivity ;
 456 #else
 457           _polyhedronNodal->setPolyhedronIndex(PolyhedronIndex);
 458           _polyhedronNodal->setFacesIndex(PolyhedronFacesIndex);
 459           _polyhedronNodal->setNodes(PolyhedronConnectivity);
 460 #endif
 461         }
 462       else
 463         {
 464           if (_typeConnectivity != MED_DESCENDING)
 465             _typeConnectivity = MED_DESCENDING; //by default it is set to MED_NODAL
 466           if (_polyhedronDescending != (MEDSKYLINEARRAY*) NULL)
 467             delete _polyhedronDescending;
 468           _polyhedronDescending = new MEDSKYLINEARRAY(NumberOfPolyhedron,ConnectivitySize,PolyhedronIndex,PolyhedronConnectivity);
 469         }
 470     }
 471   else
 472     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : _entity must be MED_CELL to set polyhedron !"));
 473 }
 474
 475 bool CONNECTIVITY::existConnectivityWithPoly (MED_EN::medConnectivity connectivityType,
 476                                               MED_EN::medEntityMesh Entity) const
 477 {
 478    if (_entity == Entity) {
 479     if ((connectivityType == MED_EN::MED_NODAL) &&
 480         (_nodal != (MEDSKYLINEARRAY*)NULL || _polygonsNodal || _polyhedronNodal))
 481       return true;
 482     if ((connectivityType == MED_EN::MED_DESCENDING) && (_descending != (MEDSKYLINEARRAY*)NULL))
 483       return true;
 484    } else if (_constituent != NULL)
 485      return _constituent->existConnectivityWithPoly(connectivityType, Entity);
 486    return false;
 487 }
 488
 489 /*! A DOCUMENTER */
 490 //------------------------------------------------------------------------------------------//
 491 void CONNECTIVITY::calculateConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity)
 492 //------------------------------------------------------------------------------------------//
 493 {
 494   MESSAGE_MED("CONNECTIVITY::calculateConnectivity");
 495
 496   // a temporary limitation due to calculteDescendingConnectivity function !!!
 497   if ((_entityDimension==3) & (Entity==MED_EDGE))
 498     throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build EDGE !");
 499
 500   if (Entity==_entity)
 501     if (ConnectivityType==MED_NODAL)
 502       calculateNodalConnectivity();
 503     else
 504       if (Entity==MED_CELL)
 505         calculateDescendingConnectivity();
 506       else
 507         throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build DESCENDING !");
 508   if (Entity!=_entity) {
 509     calculateDescendingConnectivity();
 510     if (_entityDimension == 2 || _entityDimension == 3)
 511       _constituent->calculateConnectivity(ConnectivityType,Entity);
 512   }
 513 }
 514
 515 /*!  Give, in full or no interlace mode (for nodal connectivity),
 516      descending or nodal connectivity.
 517
 518       You must give a %medEntityMesh (ie:MED_EDGE)
 519       and a %medGeometryElement (ie:MED_SEG3).  */
 520
 521 //------------------------------------------------------------//
 522 void CONNECTIVITY::updateFamily(const vector<FAMILY*>& myFamilies)
 523 //------------------------------------------------------------//
 524 {
 525   const char * LOC = "CONNECTIVITY::updateFamily(vector<FAMILY*>) ";
 526   int numberOfFamilies = myFamilies.size();
 527   if (numberOfFamilies == 0 || _constituent == NULL)
 528     return;
 529   // does we do an update ?
 530   if ((_constituent != NULL) && (_descending != NULL))
 531     return;
 532   if (myFamilies[0]->getEntity() != _constituent->getEntity())
 533     return;
 534   CONNECTIVITY * oldConstituent = _constituent;
 535   _constituent = (CONNECTIVITY *)NULL;
 536   if (oldConstituent->_nodal==NULL && oldConstituent->_polygonsNodal==NULL)
 537     throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
 538
 539   //Loc vars defined to treat polygons exactly the same as classic types. Not nice but necessary.
 540   int oldNumberOfFaceTab[2];
 541   const int * oldConstituentValueTab[2];
 542   const int * oldConstituentIndexTab[2];
 543   int * renumberingFromOldToNewTab[2];//Final mapping array between old numbers and new numbers.;
 544   int startNbOfTurnInGlobalLoop=0;
 545
 546   if (oldConstituent->_nodal != NULL)
 547   {
 548     int oldNumberOfFace = oldConstituent->_nodal->getNumberOf();
 549     oldNumberOfFaceTab[0] = oldNumberOfFace;
 550     const int * oldConstituentValue = oldConstituent->_nodal->getValue();
 551     oldConstituentValueTab[0] = oldConstituentValue;
 552     const int * oldConstituentIndex = oldConstituent->_nodal->getIndex();
 553     oldConstituentIndexTab[0] = oldConstituentIndex;
 554     int * renumberingFromOldToNew = new int [oldNumberOfFace+1];
 555     renumberingFromOldToNewTab[0] = renumberingFromOldToNew;
 556   }
 557   else //Polyg/PolyH only
 558   {
 559     renumberingFromOldToNewTab[0]=0;
 560     oldNumberOfFaceTab[0]=0;
 561     startNbOfTurnInGlobalLoop++;
 562   }
 563
 564   int oldNumberOfFacePoly = oldConstituent->getNumberOfPolygons();
 565   const int * oldConstituentValuePoly=0;
 566   const int * oldConstituentIndexPoly=0;
 567   int * renumberingFromOldToNewPoly=0;
 568
 569   int nbOfTurnInGlobalLoop=1;//Defined to know if a second search on polygons is needed.
 570   if(oldNumberOfFacePoly>0)
 571     {
 572       oldNumberOfFaceTab[1] = oldNumberOfFacePoly;
 573       oldConstituentValuePoly = oldConstituent->_polygonsNodal->getValue();
 574       oldConstituentValueTab[1] = oldConstituentValuePoly;
 575       oldConstituentIndexPoly = oldConstituent->_polygonsNodal->getIndex();
 576       oldConstituentIndexTab[1] = oldConstituentIndexPoly;
 577       renumberingFromOldToNewPoly = new int[oldNumberOfFacePoly+1];
 578       renumberingFromOldToNewTab[1] = renumberingFromOldToNewPoly;
 579       nbOfTurnInGlobalLoop++;
 580     }
 581
 582   calculateDescendingConnectivity();//perform calculation of descending connectivity to update all connectivities after having taken into account of direction of faces...
 583   _constituent->calculateReverseNodalConnectivity();//getting d-1 nodal connectivity to find get new face numbers from nodes numbers...
 584
 585   const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity(); //Common to polygons and classic geometric types
 586   const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex(); //Common to polygons and classic geometric types
 587
 588   for (int loop=startNbOfTurnInGlobalLoop; loop<nbOfTurnInGlobalLoop; loop++)
 589     {
 590       int oldNumberOfFaceLoop=oldNumberOfFaceTab[loop];
 591       const int * oldConstituentValueLoop=oldConstituentValueTab[loop];
 592       const int * oldConstituentIndexLoop= oldConstituentIndexTab[loop];
 593       int * renumberingFromOldToNewLoop=renumberingFromOldToNewTab[loop];
 594       CELLMODEL * aCELLMODEL = 0;
 595       if ( loop == 0 ) aCELLMODEL = & oldConstituent->_type[0];
 596       for(int iOldFace=0;iOldFace<oldNumberOfFaceLoop;iOldFace++)
 597         {
 598           const int *nodesOfCurrentFaceOld=oldConstituentValueLoop+oldConstituentIndexLoop[iOldFace]-1;
 599           int nbOfNodesOfCurrentFaceOld=oldConstituentIndexLoop[iOldFace+1]-oldConstituentIndexLoop[iOldFace];
 600
 601           //retrieving new number of polygon face...
 602           int sizeOfNewFaceNb=reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]]-reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]-1];
 603           int *newFaceNb=new int[sizeOfNewFaceNb];
 604           memcpy(newFaceNb,reverseFaceNodal+reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]-1]-1,sizeOfNewFaceNb*sizeof(int));//all faces candidates are in newFaceNb !!!
 605           for(int curNode=1;curNode<nbOfNodesOfCurrentFaceOld && sizeOfNewFaceNb>1;curNode++)
 606             {
 607               const int *newFacesNbForCurNode=reverseFaceNodal+reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]-1]-1;
 608               int sizeOfNewFacesNbForCurNode=reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]]-reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]-1];
 609               for(int i=0;i<sizeOfNewFaceNb;)
 610                 {
 611                   bool found=false;
 612                   for(int j=0;j<sizeOfNewFacesNbForCurNode && !found;j++)
 613                     {
 614                       if(newFacesNbForCurNode[j]==newFaceNb[i])
 615                         found=true;
 616                     }
 617                   if(found)
 618                     i++;
 619                   else
 620                     newFaceNb[i]=newFaceNb[--sizeOfNewFaceNb];
 621                 }
 622             }
 623           if(sizeOfNewFaceNb!=1)
 624             throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"unexisting face specified!!!"));
 625           renumberingFromOldToNewLoop[iOldFace]=newFaceNb[0];
 626           delete [] newFaceNb;
 627           //end of retrieving new number of polygon face...
 628           int nbOfNodesOfCurrentFaceNew;
 629           const int *nodesOfCurrentFaceNew=_constituent->getConnectivityOfAnElementWithPoly(MED_NODAL,_constituent->getEntity(),
 630                                                                                             renumberingFromOldToNewLoop[iOldFace],nbOfNodesOfCurrentFaceNew);
 631           // compare nodes of a new face and those of an old one;
 632           // for the second order elements, only vertex nodes are compared
 633           int nbOfVertices = nbOfNodesOfCurrentFaceOld;
 634           if ( aCELLMODEL ) {
 635             if ( aCELLMODEL->getNumberOfNodes() != nbOfNodesOfCurrentFaceOld ) {
 636               // type changed, find a corresponding CELLMODEL
 637               int iType = 2; // 1-st type is already used at loop beginning
 638               //while ( iOldFace + 1 >= oldConstituent->_count[ 1 + iType ]) // check next type
 639               while ( iOldFace + 1 >= oldConstituent->_count[ iType ]) // check next type
 640                 ++iType;
 641               aCELLMODEL = & oldConstituent->_type[ iType - 1 ];
 642             }
 643             nbOfVertices = aCELLMODEL->getNumberOfVertexes();
 644           }
 645           MEDMODULUSARRAY modulusArrayOld(nbOfVertices,nbOfNodesOfCurrentFaceOld,nodesOfCurrentFaceOld);
 646           int nbOfVerticesNew = nbOfVertices;
 647           if (nbOfVerticesNew > nbOfNodesOfCurrentFaceNew) nbOfVerticesNew = nbOfNodesOfCurrentFaceNew;
 648           MEDMODULUSARRAY modulusArrayNew(nbOfVerticesNew,nbOfNodesOfCurrentFaceNew,nodesOfCurrentFaceNew);
 649           int retCompareNewOld=modulusArrayNew.compare(modulusArrayOld);
 650           if(retCompareNewOld==0)
 651             throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"Uncompatible given user face with calculated existing faces"));
 652           if(retCompareNewOld==-1)
 653             invertConnectivityForAFace(renumberingFromOldToNewLoop[iOldFace],nodesOfCurrentFaceOld,loop==1);
 654         }
 655     }
 656   // Updating the Family
 657   for(vector<FAMILY*>::const_iterator iter=myFamilies.begin();iter!=myFamilies.end();iter++)
 658     (*iter)->changeElementsNbs(_constituent->getEntity(), renumberingFromOldToNewTab[0],
 659                                oldNumberOfFaceTab[0], renumberingFromOldToNewTab[1]);
 660
 661   // FIX PAL13414:
 662   if ( _constituent && !_constituent->_constituent ) {
 663     _constituent->_constituent = oldConstituent->_constituent;
 664     oldConstituent->_constituent = NULL;
 665   }
 666   // END FIX PAL13414:
 667
 668   delete oldConstituent;
 669   delete [] renumberingFromOldToNewTab[0];
 670   if (oldNumberOfFacePoly > 0)
 671     delete [] renumberingFromOldToNewPoly;
 672   return;
 673 }
 674
 675 //------------------------------------------------------------------------------------------------------------------//
 676 const int * MEDMEM::CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type)
 677 //------------------------------------------------------------------------------------------------------------------//
 678 {
 679   const char * LOC = "CONNECTIVITY::getConnectivity";
 680   //  BEGIN_OF_MED(LOC);
 681
 682   MEDSKYLINEARRAY * Connectivity;
 683   if (Entity==_entity) {
 684
 685     if (ConnectivityType==MED_NODAL)
 686       {
 687         calculateNodalConnectivity();
 688         Connectivity=_nodal;
 689       }
 690     else
 691       {
 692         calculateDescendingConnectivity();
 693         Connectivity=_descending;
 694       }
 695
 696     if (Connectivity!=NULL)
 697       if (Type==MED_ALL_ELEMENTS)
 698         return Connectivity->getValue();
 699       else {
 700         for (int i=0; i<_numberOfTypes; i++)
 701           if (_geometricTypes[i]==Type)
 702             //return Connectivity->getI(i+1);
 703             return Connectivity->getI(_count[i]);
 704         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
 705       }
 706     else
 707       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
 708   } else
 709     if (_constituent != NULL)
 710       return _constituent->getConnectivity(ConnectivityType,Entity,Type);
 711
 712   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 713 }
 714
 715 //------------------------------------------------------------------------------------------------------------------//
 716 int CONNECTIVITY::getConnectivityLength(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type)
 717 //------------------------------------------------------------------------------------------------------------------//
 718 {
 719   const char * LOC = "CONNECTIVITY::getConnectivity";
 720   BEGIN_OF_MED(LOC);
 721
 722   MEDSKYLINEARRAY * Connectivity;
 723   if (Entity==_entity) {
 724
 725     if (ConnectivityType==MED_NODAL)
 726       {
 727         calculateNodalConnectivity();
 728         Connectivity=_nodal;
 729       }
 730     else
 731       {
 732         calculateDescendingConnectivity();
 733         Connectivity=_descending;
 734       }
 735
 736     if (Connectivity!=NULL)
 737       if (Type==MED_ALL_ELEMENTS)
 738         return Connectivity->getLength();
 739       else {
 740         for (int i=0; i<_numberOfTypes; i++)
 741           if (_geometricTypes[i]==Type)
 742             // issue 19983
 743             //return (_count[i+1]-_count[i])*getType(Type).getNumberOfNodes();
 744           {
 745             const int *ind=Connectivity->getIndex();
 746             return ind[_count[i+1]-1]-ind[_count[i]-1];
 747           }
 748         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
 749       }
 750     else
 751       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
 752   }
 753   else
 754     if (_constituent != NULL)
 755       return _constituent->getConnectivityLength(ConnectivityType,Entity,Type);
 756
 757   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 758 }
 759
 760 /*!  Give morse index array to use with
 761      getConnectivity(MED_FULL_INTERLACE,mode,entity,MED_ALL_ELEMENTS).
 762
 763       Each value give start index for corresponding entity in connectivity array.
 764
 765       Example : i-th element, j-th node of it :
 766       - In C mode : Connectivity[ConnectivityIndex[i]-1+j-1]
 767       - In fortran mode : Connectivity[ConnectivityIndex[i]+j] */
 768 //-----------------------------------------------------------------------------------------------//
 769 const int * CONNECTIVITY::getConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity)
 770 //-----------------------------------------------------------------------------------------------//
 771 {
 772   const char * LOC = "CONNECTIVITY::getConnectivityIndex";
 773
 774   MEDSKYLINEARRAY * Connectivity;
 775   if (Entity==_entity) {
 776
 777     if (ConnectivityType==MED_NODAL)
 778       Connectivity=_nodal;
 779     else
 780       Connectivity=_descending;
 781
 782     if (Connectivity!=NULL)
 783       return Connectivity->getIndex();
 784     else
 785       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
 786
 787   } else
 788     if (_constituent != NULL)
 789       return _constituent->getConnectivityIndex(ConnectivityType,Entity);
 790
 791   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 792 }
 793
 794
 795 //-------------------------------------------------------------//
 796 const int* CONNECTIVITY::getPolygonsConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity)
 797 //-------------------------------------------------------------//
 798 {
 799   const char* LOC = "CONNECTIVITY::getPolygonsConnectivity";
 800   BEGIN_OF_MED(LOC);
 801
 802   MEDSKYLINEARRAY* Connectivity;
 803   if (Entity == _entity)
 804     {
 805       if (ConnectivityType == MED_NODAL)
 806         {
 807           calculateNodalConnectivity();
 808           Connectivity = _polygonsNodal;
 809         }
 810       else
 811         {
 812           calculateDescendingConnectivity();
 813           Connectivity = _polygonsDescending;
 814         }
 815       if (Connectivity != NULL)
 816         return Connectivity->getValue();
 817       else
 818         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polygons Connectivity not defined !"));
 819     }
 820   else
 821     if (_constituent != NULL)
 822       return _constituent->getPolygonsConnectivity(ConnectivityType, Entity);
 823   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 824 }
 825
 826
 827 //-------------------------------------------------------------//
 828 const int* CONNECTIVITY::getPolygonsConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity)
 829 //-------------------------------------------------------------//
 830 {
 831   const char* LOC = "CONNECTIVITY::getPolygonsConnectivityIndex";
 832   BEGIN_OF_MED(LOC);
 833
 834   MEDSKYLINEARRAY* Connectivity;
 835   if (Entity == _entity)
 836     {
 837       if (ConnectivityType == MED_NODAL)
 838         {
 839           //      calculateNodalConnectivity();
 840           Connectivity = _polygonsNodal;
 841         }
 842       else
 843         {
 844           //      calculateDescendingConnectivity();
 845           Connectivity = _polygonsDescending;
 846         }
 847       if (Connectivity != NULL)
 848         return Connectivity->getIndex();
 849       else
 850         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polygons Connectivity not defined !"));
 851     }
 852   else
 853     if (_constituent != NULL)
 854       return _constituent->getPolygonsConnectivityIndex(ConnectivityType, Entity);
 855   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
 856 }
 857
 858
 859 /*! We suppose in this method that nodal and descending connectivities
 860   are coherent.*/
 861 //-------------------------------------------------------------//
 862 int CONNECTIVITY::getNumberOfPolygons(MED_EN::medEntityMesh Entity) const
 863 //-------------------------------------------------------------//
 864 {
 865   if(Entity==MED_ALL_ENTITIES || Entity==_entity )
 866     {
 867       if (_polygonsNodal != (MEDSKYLINEARRAY*) NULL)
 868         return _polygonsNodal->getNumberOf();
 869       else if (_polygonsDescending != (MEDSKYLINEARRAY*) NULL)
 870         return _polygonsDescending->getNumberOf();
 871       else
 872         return 0;
 873     }
 874   else
 875     {
 876       if (_constituent == (CONNECTIVITY*) NULL)
 877         throw MEDEXCEPTION("getNumberOfPolygons : Entity not found !");
 878       return _constituent->getNumberOfPolygons(Entity);
 879     }
 880 }
 881
 882
 883 //--------------------------------------------------------------//
 884 const int* CONNECTIVITY::getPolyhedronConnectivity(medConnectivity ConnectivityType) const
 885 //--------------------------------------------------------------//
 886 {
 887   const char* LOC = "CONNECTIVITY::getPolyhedronConnectivity";
 888   BEGIN_OF_MED(LOC);
 889
 890   if (_entity == MED_CELL) //polyhedron can only be MED_CELL
 891     {
 892       if (ConnectivityType == MED_NODAL)
 893         {
 894           ((CONNECTIVITY *)(this))->calculateNodalConnectivity();
 895           if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL) {
 896 //CCRT
 897 #if defined(IRIX64) || defined(OSF1) || defined(VPP5000) || defined(PCLINUX64)
 898             int i ;
 899             const MED_EN::med_int * tmp_PolyhedronConnectivity = _polyhedronNodal->getNodes();
 900             int * PolyhedronConnectivity = new int[_polyhedronNodal->getNumberOfNodes()] ;
 901             for ( i = 0 ; i < _polyhedronNodal->getNumberOfNodes() ; i++ )
 902                PolyhedronConnectivity[i] = tmp_PolyhedronConnectivity[i] ;
 903 //CCRT : return of a copy of PolyhedronConnectivity
 904             return PolyhedronConnectivity ;
 905 #else
 906             return _polyhedronNodal->getNodes();
 907 #endif
 908           }
 909           else
 910             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polyhedron Nodal Connectivity not defined !"));
 911         }
 912       else
 913         {
 914           ((CONNECTIVITY *)(this))->calculateDescendingConnectivity();
 915           if (_polyhedronDescending != (MEDSKYLINEARRAY*) NULL) {
 916             return _polyhedronDescending->getValue();
 917           }
 918           else
 919             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polyhedron Descending Connectivity not defined !"));
 920         }
 921     }
 922   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : No polyhedron for _entity different from MED_CELL !"));
 923 }
 924
 925
 926 //---------------------------------------------------------------//
 927 const int* CONNECTIVITY::getPolyhedronFacesIndex() const
 928 //---------------------------------------------------------------//
 929 {
 930   const char* LOC = "CONNECTIVITY::getPolyhedronFacesIndex";
 931   BEGIN_OF_MED(LOC);
 932
 933   if (_entity == MED_CELL) //polyhedron can only be MED_CELL
 934     {
 935       //      calculateNodalConnectivity();
 936       if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL) {
 937 //CCRT
 938 #if defined(IRIX64) || defined(OSF1) || defined(VPP5000) || defined(PCLINUX64)
 939             int i ;
 940             const MED_EN::med_int * tmp_PolyhedronFacesIndex = _polyhedronNodal->getFacesIndex();
 941             int * PolyhedronFacesIndex = new int[_polyhedronNodal->getNumberOfFaces()+1] ;
 942             for ( i = 0 ; i < _polyhedronNodal->getNumberOfFaces()+1 ; i++ )
 943                PolyhedronFacesIndex[i] = tmp_PolyhedronFacesIndex[i] ;
 944 //CCRT : return of a copy of PolyhedronFacesIndex
 945             return PolyhedronFacesIndex ;
 946 #else
 947         return _polyhedronNodal->getFacesIndex();
 948 #endif
 949       }
 950       else
 951         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : No Polyhedron in that Connectivity !"));
 952     }
 953   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : No polyhedron for _entity different from MED_CELL !"));
 954 }
 955
 956
 957 //--------------------------------------------------------------//
 958 const int* CONNECTIVITY::getPolyhedronIndex(medConnectivity ConnectivityType) const
 959 //--------------------------------------------------------------//
 960 {
 961   const char* LOC = "CONNECTIVITY::getPolyhedronIndex";
 962   BEGIN_OF_MED(LOC);
 963
 964   if (_entity == MED_CELL) //polyhedron can only be MED_CELL
 965     {
 966       if (ConnectivityType == MED_NODAL)
 967         {
 968           //      calculateNodalConnectivity();
 969           if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL) {
 970 //CCRT
 971 #if defined(IRIX64) || defined(OSF1) || defined(VPP5000) || defined(PCLINUX64)
 972             int i ;
 973             const MED_EN::med_int * tmp_PolyhedronIndex = _polyhedronNodal->getPolyhedronIndex();
 974             int * PolyhedronIndex = new int[_polyhedronNodal->getNumberOfPolyhedron()+1] ;
 975             for ( i = 0 ; i < _polyhedronNodal->getNumberOfPolyhedron()+1 ; i++ )
 976                PolyhedronIndex[i] = tmp_PolyhedronIndex[i] ;
 977 //CCRT : return of a copy of PolyhedronIndex
 978             return PolyhedronIndex ;
 979 #else
 980             return _polyhedronNodal->getPolyhedronIndex();
 981 #endif
 982           }
 983           else
 984             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polyhedron Nodal Connectivity not defined !"));
 985         }
 986       else
 987         {
 988           //      calculateDescendingConnectivity();
 989           if (_polyhedronDescending != (MEDSKYLINEARRAY*) NULL)
 990             return _polyhedronDescending->getIndex();
 991           else
 992             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Polyhedron Descending Connectivity not defined !"));
 993         }
 994     }
 995   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : No polyhedron for _entity different from MED_CELL !"));
 996 }
 997
 998
 999 /*! We suppose in this method that nodal and descending connectivities
1000   are coherent.*/
1001 //-------------------------------------------------------------//
1002 int CONNECTIVITY::getNumberOfPolyhedronFaces() const
1003 //-------------------------------------------------------------//
1004 {
1005   //  if (_polyhedronNodal == (POLYHEDRONARRAY*) NULL)
1006   //    calculateNodalConnectivity();
1007   if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
1008     return _polyhedronNodal->getNumberOfFaces();
1009   else
1010     return 0;
1011 }
1012
1013
1014 /*! We suppose in this method that nodal and descending connectivities
1015   are coherent.*/
1016 //--------------------------------------------------------------//
1017 int CONNECTIVITY::getNumberOfPolyhedron() const
1018 //--------------------------------------------------------------//
1019 {
1020   if (_polyhedronNodal != (POLYHEDRONARRAY*) NULL)
1021     return _polyhedronNodal->getNumberOfPolyhedron();
1022   else if (_polyhedronDescending != (MEDSKYLINEARRAY*) NULL)
1023     return _polyhedronDescending->getNumberOf();
1024   else
1025     return 0;
1026 }
1027
1028
1029 /*! A DOCUMENTER */
1030 //--------------------------------------------------------------//
1031 const CELLMODEL & CONNECTIVITY::getType(medGeometryElement Type) const
1032 //--------------------------------------------------------------//
1033 {
1034
1035   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
1036     throw MEDEXCEPTION("CONNECTIVITY::getType : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE !");
1037   for (int i=0; i<_numberOfTypes; i++)
1038     if (_geometricTypes[i]==Type)
1039       return _type[i];
1040   throw MEDEXCEPTION("CONNECTIVITY::getType :  medGeometryElement not found !");
1041 }
1042
1043 /*!  Returns the number of elements of type %medGeometryElement.
1044      Note : not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
1045 //------------------------------------------------------------------------//
1046 int CONNECTIVITY::getNumberOfNodesInType(medGeometryElement Type) const
1047 //------------------------------------------------------------------------//
1048 {
1049   const char * LOC = "CONNECTIVITY::getNumberOfNodesInType";
1050   BEGIN_OF_MED(LOC);
1051
1052   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
1053     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!"));
1054   for (int i=0; i<_numberOfTypes; i++)
1055     if (_geometricTypes[i]==Type)
1056       return _type[i].getNumberOfNodes();
1057   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
1058 }
1059
1060 /*!  Returns the number of geometric sub cells of %medGeometryElement type.
1061 not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
1062 //------------------------------------------------------------------------//
1063 int CONNECTIVITY::getNumberOfSubCellInType(medGeometryElement Type) const
1064 //------------------------------------------------------------------------//
1065 {
1066   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
1067     throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!");
1068   for (int i=0; i<_numberOfTypes; i++)
1069     if (_geometricTypes[i]==Type)
1070       return _type[i].getNumberOfConstituents(1);
1071   throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement not found !");
1072 }
1073
1074 /*!  Returns the number of elements of type %medGeometryElement.
1075
1076      Note :
1077       - Implemented for MED_ALL_ELEMENTS
1078       - Not implemented for MED_ALL_ENTITIES (A VERIFIER)
1079       - Not implemented for MED_NONE */
1080 //-----------------------------------------------------------------------------------//
1081 int CONNECTIVITY::getNumberOf(medEntityMesh Entity, medGeometryElement Type) const
1082 //-----------------------------------------------------------------------------------//
1083 {
1084   if (Entity==_entity) {
1085     if (Type==MED_EN::MED_NONE)
1086       return 0; // not defined !
1087     //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_NONE"));
1088     if (!existConnectivity(MED_NODAL,Entity) && !existConnectivity(MED_DESCENDING,Entity)) return 0; //case with only polygons for example
1089     if (Type==MED_EN::MED_ALL_ELEMENTS)
1090       return _count[_numberOfTypes]-1;
1091     for (int i=0; i<_numberOfTypes; i++)
1092       if (_geometricTypes[i]==Type)
1093         return (_count[i+1] - _count[i]);
1094     //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
1095   } else
1096     if (_constituent != NULL)
1097       return _constituent->getNumberOf(Entity,Type);
1098
1099   return 0; // valid if they are nothing else !
1100   //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
1101 }
1102
1103 /*! A DOCUMENTER */
1104 //--------------------------------------------------------------//
1105 const int* CONNECTIVITY::getValue(medConnectivity TypeConnectivity,
1106                                 medGeometryElement Type)
1107 //--------------------------------------------------------------//
1108 {
1109   if (TypeConnectivity == MED_NODAL)
1110     {
1111           calculateNodalConnectivity();
1112           if (Type==MED_ALL_ELEMENTS)
1113             return _nodal->getValue();
1114           for (int i=0; i<_numberOfTypes; i++)
1115             if (_geometricTypes[i]==Type)
1116               return _nodal->getI(_count[i]);
1117     }
1118   else
1119     {
1120       calculateDescendingConnectivity();
1121       if (Type==MED_ALL_ELEMENTS)
1122         return _descending->getValue();
1123       for (int i=0; i<_numberOfTypes; i++)
1124         if (_geometricTypes[i]==Type)
1125           return _descending->getI(_count[i]);
1126     }
1127   throw MEDEXCEPTION("Not found");
1128 }
1129
1130 /*! A DOCUMENTER */
1131 //---------------------------------------------------------------------//
1132 const int* CONNECTIVITY:: getValueIndex(medConnectivity TypeConnectivity)
1133 //---------------------------------------------------------------------//
1134 {
1135   if (TypeConnectivity == MED_NODAL)
1136     {
1137       calculateNodalConnectivity();
1138       return _nodal->getIndex();
1139     }
1140   else
1141     {
1142       calculateDescendingConnectivity();
1143       return _descending->getIndex();
1144     }
1145 }
1146
1147 /*! Not yet implemented */
1148 //----------------------------------------------//
1149 const int* CONNECTIVITY:: getNeighbourhood() const
1150 //----------------------------------------------//
1151 {
1152   throw MEDEXCEPTION("CONNECTIVITY::getNeighbourhood : Not implemented");
1153 }
1154
1155 /*! Returns an array which contains, for each node, all cells
1156     arround it. */
1157 //-------------------------------------------------//
1158 const int* CONNECTIVITY::getReverseNodalConnectivity()
1159 //-------------------------------------------------//
1160 {
1161   calculateReverseNodalConnectivity();
1162   return _reverseNodalConnectivity->getValue();
1163 }
1164
1165 /*!  Give index array to use with getReverseConnectivity(MED_NODAL).
1166      It is unusefull with MED_DESCENDING mode, because we have allways two cells.  */
1167 //-------------------------------------------------------//
1168 const int* CONNECTIVITY::getReverseNodalConnectivityIndex()
1169 //-------------------------------------------------------//
1170 {
1171   calculateReverseNodalConnectivity();
1172   return _reverseNodalConnectivity->getIndex();
1173 }
1174
1175 /*! Returns an array which contains, for each face (or edge),
1176     the 2 cells of each side. First is cell which face normal is outgoing.
1177     arround it. */
1178 //------------------------------------------------------//
1179 const int* CONNECTIVITY::getReverseDescendingConnectivity()
1180 //------------------------------------------------------//
1181 {
1182   // it is in _constituent connectivity only if we are in MED_CELL
1183   // (we could not for instance calculate face-edge connectivity !)
1184   if (_entity!=MED_CELL)
1185     throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivity : Error Only in MED_CELL connectivity");
1186
1187   // we need descending connectivity
1188   calculateDescendingConnectivity();
1189         if (_reverseDescendingConnectivity==NULL)
1190                 _reverseDescendingConnectivity=_descending->makeReverseArray();
1191
1192   return _reverseDescendingConnectivity->getValue();
1193 }
1194
1195 /*! calculate the reverse descending Connectivity
1196     and returns the index  ( A DOCUMENTER MIEUX)*/
1197 //-----------------------------------------------------------//
1198 const int* CONNECTIVITY::getReverseDescendingConnectivityIndex()
1199 //-----------------------------------------------------------//
1200 {
1201   // it is in _constituent connectivity only if we are in MED_CELL
1202   if (_entity!=MED_CELL)
1203     throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivityIndex : Error Only in MED_CELL connectivity");
1204
1205   // we need descending connectivity
1206   calculateDescendingConnectivity();
1207   return _reverseDescendingConnectivity->getIndex();
1208 }
1209
1210 /*! A DOCUMENTER (et a finir ???) */
1211 //--------------------------------------------//
1212 void CONNECTIVITY::calculateNodalConnectivity()
1213 //--------------------------------------------//
1214 {
1215   if (_nodal==NULL && _polygonsNodal==NULL && _polyhedronNodal==NULL)
1216     {
1217       if (_descending==NULL && _polygonsDescending==NULL && _polyhedronDescending==NULL)
1218         throw MEDEXCEPTION("CONNECTIVITY::calculateNodalConnectivity : No connectivity found !");
1219       // calculate _nodal, _polygonsNodal and _polyhedronNodal from _descending, _polygonsDescending and _polyhedronDescending
1220     }
1221 }
1222
1223
1224
1225 /*! If not yet done, calculate the nodal Connectivity
1226     and the reverse nodal Connectivity*/
1227 //---------------------------------------------------//
1228 void CONNECTIVITY::calculateReverseNodalConnectivity()
1229 //---------------------------------------------------//
1230 {
1231   const char* LOC = "CONNECTIVITY::calculateReverseNodalConnectivity : ";
1232   BEGIN_OF_MED(LOC);
1233
1234   SCRUTE_MED(_nodal);
1235   SCRUTE_MED(_reverseNodalConnectivity);
1236
1237   if (_nodal==NULL)
1238     calculateNodalConnectivity();
1239
1240   if(_reverseNodalConnectivity==NULL) {
1241
1242     int node_number = 0;
1243     vector <vector <int> > reverse_connectivity;
1244     reverse_connectivity.resize(_numberOfNodes+1);
1245
1246   // Treat all cells types
1247
1248     for (int j = 0; j < _numberOfTypes; j++)
1249       {
1250         // node number of the cell type
1251         node_number = _type[j].getNumberOfNodes();
1252         // treat all cells of a particular type
1253         for (int k = _count[j]; k < _count[j+1]; k++)
1254           // treat all nodes of the cell type
1255           for (int local_node_number = 1; local_node_number < node_number+1; local_node_number++)
1256             {
1257               int global_node = _nodal->getIJ(k,local_node_number);
1258               reverse_connectivity[global_node].push_back(k);
1259             }
1260       }
1261
1262     if(_polygonsNodal != NULL && _entityDimension==2)
1263       {
1264         int nbOfPolygons=_polygonsNodal->getNumberOf();
1265         int offset=getNumberOf(_entity,MED_ALL_ELEMENTS);
1266         const int *index=_polygonsNodal->getIndex();
1267         const int *value=_polygonsNodal->getValue();
1268         for (int local_polyg_number = 0; local_polyg_number < nbOfPolygons; local_polyg_number++)
1269           {
1270             for(int i=index[local_polyg_number];i<index[local_polyg_number+1];i++)
1271               reverse_connectivity[value[i-1]].push_back(offset+local_polyg_number+1);
1272           }
1273       }
1274
1275     if(_polyhedronNodal != NULL && _entityDimension==3)
1276       {
1277         int nbOfPolyhedra=_polyhedronNodal->getNumberOfPolyhedron();
1278         int offset=getNumberOf(_entity,MED_ALL_ELEMENTS);
1279         for (int local_polyh_number = 0; local_polyh_number < nbOfPolyhedra; local_polyh_number++)
1280             {
1281               int nbOfNodes;
1282               int global_polyh_number=offset+local_polyh_number+1;
1283               int *nodes=getNodesOfPolyhedron(global_polyh_number,nbOfNodes);
1284               for(int i=0;i<nbOfNodes;i++)
1285                 reverse_connectivity[nodes[i]].push_back(global_polyh_number);
1286               delete [] nodes;
1287             }
1288       }
1289
1290     // Full reverse_nodal_connectivity and reverse_nodal_connectivity_index
1291
1292   //calculate size of reverse_nodal_connectivity array
1293     int size_reverse_nodal_connectivity  = 0;
1294     for (int i = 1; i < _numberOfNodes+1; i++)
1295       size_reverse_nodal_connectivity += reverse_connectivity[i].size();
1296
1297     int * reverse_nodal_connectivity_index = new int[_numberOfNodes+1];
1298     int * reverse_nodal_connectivity = new int[size_reverse_nodal_connectivity];
1299
1300     reverse_nodal_connectivity_index[0] = 1;
1301     for (int i = 1; i < _numberOfNodes+1; i++)
1302       {
1303         int size = reverse_connectivity[i].size();
1304         reverse_nodal_connectivity_index[i] = reverse_nodal_connectivity_index[i-1] + size;
1305         for (int j = 0; j < size; j++)
1306           reverse_nodal_connectivity[reverse_nodal_connectivity_index[i-1]-1+j]= reverse_connectivity[i][j];
1307       }
1308
1309     _reverseNodalConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity,
1310                                                     reverse_nodal_connectivity_index,
1311                                                     reverse_nodal_connectivity,true);
1312   }
1313   END_OF_MED(LOC);
1314 }
1315
1316 /*! If not yet done, calculate the Descending Connectivity */
1317 //-------------------------------------------------//
1318 void CONNECTIVITY::calculateDescendingConnectivity()
1319   //-------------------------------------------------//
1320 {
1321   const char * LOC = "CONNECTIVITY::calculateDescendingConnectivity() : ";
1322   BEGIN_OF_MED(LOC);
1323
1324   if (_descending==NULL && _polygonsDescending==NULL && _polyhedronDescending==NULL)
1325   {
1326     if (_nodal==NULL && _polygonsNodal==NULL && _polyhedronNodal==NULL)
1327     {
1328       MESSAGE_MED(LOC<<"No connectivity found !");
1329       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
1330     }
1331     // calcul _descending from _nodal
1332     // we need CONNECTIVITY for constituent
1333
1334     if (_constituent != NULL && _constituent->_nodal != NULL)
1335     {
1336       calculatePartialDescendingConnectivity();
1337     }
1338     else
1339     {
1340       calculateFullDescendingConnectivity(_entity);
1341     }
1342   }
1343 }
1344
1345 /*! If not yet done, calculate the full Descending Connectivity */
1346 //-------------------------------------------------//
1347 void CONNECTIVITY::calculateFullDescendingConnectivity(MED_EN::medEntityMesh Entity)
1348 //-------------------------------------------------//
1349 {
1350   const char * LOC = "CONNECTIVITY::calculateFullDescendingConnectivity() : ";
1351   BEGIN_OF_MED(LOC);
1352   if (_entity != Entity)
1353   {
1354     if (_constituent == NULL)
1355       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Entity not found !"));
1356     _constituent->calculateFullDescendingConnectivity(Entity);
1357   }
1358   else
1359   {
1360     if (_descending==NULL && _polygonsDescending==NULL && _polyhedronDescending==NULL ||
1361         _isDescendingConnectivityPartial )
1362     {
1363       if (_nodal==NULL && _polygonsNodal==NULL && _polyhedronNodal==NULL)
1364       {
1365         MESSAGE_MED(LOC<<"No connectivity found !");
1366         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
1367       }
1368       delete _constituent;
1369       _constituent=0;
1370       if (_entityDimension == 3)
1371         _constituent = new CONNECTIVITY(MED_FACE);
1372       else if (_entityDimension == 2)
1373         _constituent = new CONNECTIVITY(MED_EDGE);
1374       else
1375         {
1376           MESSAGE_MED(LOC<<"We are in 1D");
1377           return;
1378         }
1379
1380       _constituent->_typeConnectivity = MED_NODAL;
1381       _constituent->_numberOfNodes = _numberOfNodes;
1382       // foreach cells, we built array of constituent
1383       int DescendingSize = 0;
1384       for(int i=0; i<_numberOfTypes; i++)
1385         DescendingSize+=(_count[i+1]-_count[i])*_type[i].getNumberOfConstituents(1);
1386       //_descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,DescendingSize);
1387       //const int * descend_connectivity = _descending->getValue();
1388       int * descend_connectivity = new int[DescendingSize];
1389       for (int i=0; i<DescendingSize; i++)
1390         descend_connectivity[i]=0;
1391       //const int * descend_connectivity_index = _descending->getIndex();
1392       int * descend_connectivity_index = new int[_count[_numberOfTypes]];
1393       if(_numberOfTypes>0)
1394         descend_connectivity_index[0]=1;
1395
1396
1397       map<medGeometryElement,int> eltsCounter;
1398       medGeometryElement* ConstituentsTypes = new medGeometryElement[2];
1399       ConstituentsTypes[0]=MED_NONE;
1400       ConstituentsTypes[1]=MED_NONE;
1401       int * NumberOfConstituentsForeachType = new int[2];
1402       NumberOfConstituentsForeachType[0]=0;
1403       NumberOfConstituentsForeachType[1]=0;
1404       map<medGeometryElement,int>::iterator status;
1405       for(int i=0; i<_numberOfTypes; i++)
1406         {
1407           // initialize descend_connectivity_index array :
1408           int NumberOfConstituents = _type[i].getNumberOfConstituents(1);
1409           for (int j=_count[i];j<_count[i+1];j++)
1410             {
1411               descend_connectivity_index[j]=descend_connectivity_index[j-1]+NumberOfConstituents;
1412               // compute number of constituent of all cell for each type
1413               for(int k=1;k<NumberOfConstituents+1;k++)
1414                 {
1415                   medGeometryElement MEDType = _type[i].getConstituentType(1,k);
1416                   status=eltsCounter.find(MEDType);
1417                   if(status!=eltsCounter.end())
1418                     (*status).second++;
1419                   else
1420                     eltsCounter[MEDType]=1;
1421                 }
1422             }
1423         }
1424       if(eltsCounter.size()>2)
1425         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" Descending connectivity does not support more than 2 types."));
1426       status=eltsCounter.begin();
1427       if(!eltsCounter.empty())
1428       {
1429         ConstituentsTypes[0]=(*status).first; NumberOfConstituentsForeachType[0]=(*status).second;
1430         status++;
1431         if(status!=eltsCounter.end())
1432           { ConstituentsTypes[1]=(*status).first; NumberOfConstituentsForeachType[1]=(*status).second; }
1433       }
1434       // we could built _constituent !
1435       int TotalNumberOfConstituents = NumberOfConstituentsForeachType[0]+NumberOfConstituentsForeachType[1];
1436       int TotalNumberOfNodes = NumberOfConstituentsForeachType[0]*(ConstituentsTypes[0]%100)+NumberOfConstituentsForeachType[1]*(ConstituentsTypes[1]%100);
1437
1438       //_constituent->_nodal = new MEDSKYLINEARRAY(TotalNumberOfConstituents,TotalNumberOfNodes);
1439
1440       // we use _constituent->_nodal
1441       //const int * ConstituentNodalConnectivity = _constituent->_nodal->getValue();
1442       //const int * ConstituentNodalConnectivityIndex = _constituent->_nodal->getIndex();
1443       int * ConstituentNodalConnectivity = new int[TotalNumberOfNodes];
1444       int * ConstituentNodalConnectivityIndex = new int[TotalNumberOfConstituents+1];
1445       ConstituentNodalConnectivityIndex[0]=1;
1446
1447       _constituent->_entityDimension = _entityDimension-1;
1448       if(ConstituentsTypes[0]==MED_NONE && ConstituentsTypes[1]==MED_NONE && _numberOfTypes==0)
1449         _constituent->_numberOfTypes = 0;
1450       else if (ConstituentsTypes[1]==MED_NONE)
1451         _constituent->_numberOfTypes = 1;
1452       else
1453         _constituent->_numberOfTypes = 2;
1454       _constituent->_geometricTypes = new medGeometryElement[_constituent->_numberOfTypes];
1455       _constituent->_type = new CELLMODEL[_constituent->_numberOfTypes];
1456 //CCRT      _constituent->_count = new med_int[_constituent->_numberOfTypes+1];
1457       if(_constituent->_count)
1458         delete [] _constituent->_count;
1459       _constituent->_count = new int[_constituent->_numberOfTypes+1];
1460       _constituent->_count[0]=1;
1461       med_int* tmp_NumberOfConstituentsForeachType = new med_int[_constituent->_numberOfTypes+1];
1462       tmp_NumberOfConstituentsForeachType[0]=0; // to count constituent of each type
1463       for (int i=0; i<_constituent->_numberOfTypes;i++) {
1464         _constituent->_geometricTypes[i]=ConstituentsTypes[i];
1465         _constituent->_count[i+1]=_constituent->_count[i]+NumberOfConstituentsForeachType[i];
1466         CELLMODEL Type(ConstituentsTypes[i]);
1467         _constituent->_type[i]=Type;
1468         tmp_NumberOfConstituentsForeachType[i+1]=NumberOfConstituentsForeachType[i];
1469         for (int j=tmp_NumberOfConstituentsForeachType[i]; j<tmp_NumberOfConstituentsForeachType[i+1]+tmp_NumberOfConstituentsForeachType[i]; j++)
1470           ConstituentNodalConnectivityIndex[j+1]=ConstituentNodalConnectivityIndex[j]+(ConstituentsTypes[i]%100);
1471       }
1472       delete [] ConstituentsTypes;
1473       delete [] NumberOfConstituentsForeachType;
1474
1475       // we need reverse nodal connectivity
1476       if (! _reverseNodalConnectivity)
1477         calculateReverseNodalConnectivity();
1478       const int * ReverseNodalConnectivityIndex = _reverseNodalConnectivity->getIndex();
1479       const int * ReverseNodalConnectivityValue = _reverseNodalConnectivity->getValue();
1480
1481       // array to keep reverse descending connectivity
1482       int * ReverseDescendingConnectivityValue = new int[TotalNumberOfConstituents*2];
1483
1484       int TotalNumberOfSubCell = 0;
1485       for (int i=0; i<_numberOfTypes; i++)
1486         { // loop on all cell type
1487           CELLMODEL Type = _type[i];
1488           //      int NumberOfNodesPerCell = Type.getNumberOfNodes();
1489           int NumberOfConstituentPerCell = Type.getNumberOfConstituents(1);
1490           for (int j=_count[i]; j<_count[i+1]; j++) // we loop on all cell of this type
1491             for (int k=1; k<=NumberOfConstituentPerCell; k++)
1492               { // we loop on all sub cell of it
1493                 if (descend_connectivity[descend_connectivity_index[j-1]+k-2]==0) {
1494                   // it is a new sub cell !
1495                   //        TotalNumberOfSubCell++;
1496                   // Which type ?
1497                   if (Type.getConstituentType(1,k)==_constituent->_geometricTypes[0]){
1498                     tmp_NumberOfConstituentsForeachType[0]++;
1499                     TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[0];
1500                   } else {
1501                     tmp_NumberOfConstituentsForeachType[1]++;
1502                     TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[1];
1503                   }
1504                   //we have maximum two types
1505
1506                   descend_connectivity[descend_connectivity_index[j-1]+k-2]=TotalNumberOfSubCell;
1507                   ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2]=j;
1508                   int NumberOfNodesPerConstituent = Type.getConstituentType(1,k)%100;
1509
1510                   int * NodesLists = new int[NumberOfNodesPerConstituent];
1511                   for (int l=0; l<NumberOfNodesPerConstituent; l++) {
1512                     NodesLists[l]=_nodal->getIJ(j,Type.getNodeConstituent(1,k,l+1));
1513                     ConstituentNodalConnectivity[ConstituentNodalConnectivityIndex[TotalNumberOfSubCell-1]-1+l]=NodesLists[l];
1514                   }
1515                   // we use reverse_nodal_connectivity to find the other element which contain this sub cell
1516
1517                   // all elements which contains first node of sub cell :
1518                   int ReverseNodalConnectivityIndex_0 = ReverseNodalConnectivityIndex[NodesLists[0]-1];
1519                   int ReverseNodalConnectivityIndex_1 = getIndexOfEndClassicElementInReverseNodal(ReverseNodalConnectivityValue,ReverseNodalConnectivityIndex,NodesLists[0]-1);
1520                   //ReverseNodalConnectivityIndex[NodesLists[0]];
1521                   int NumberOfCellsInList = ReverseNodalConnectivityIndex_1-ReverseNodalConnectivityIndex_0;
1522
1523                   if (NumberOfCellsInList > 0)
1524                     { // we could have no element !
1525                       int * CellsList = new int[NumberOfCellsInList];
1526                       for (int l=ReverseNodalConnectivityIndex_0; l<ReverseNodalConnectivityIndex_1; l++)
1527                         CellsList[l-ReverseNodalConnectivityIndex_0]=ReverseNodalConnectivityValue[l-1];
1528
1529                       // foreach node in sub cell, we search elements which are in common
1530                       // at the end, we must have only one !
1531
1532                       for (int l=1; l<NumberOfNodesPerConstituent; l++) {
1533                         int NewNumberOfCellsInList = 0;
1534                         int * NewCellsList = new int[NumberOfCellsInList];
1535                         for (int l1=0; l1<NumberOfCellsInList; l1++)
1536                           for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<getIndexOfEndClassicElementInReverseNodal(ReverseNodalConnectivityValue,ReverseNodalConnectivityIndex,NodesLists[l]-1); l2++)
1537                           //for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<ReverseNodalConnectivityIndex[NodesLists[l]]; l2++)
1538                             {
1539                               if (CellsList[l1]<ReverseNodalConnectivityValue[l2-1])
1540                                 // increasing order : CellsList[l1] are not in elements list of node l
1541                                 break;
1542                               if ((CellsList[l1]==ReverseNodalConnectivityValue[l2-1])&(CellsList[l1]!=j)) {
1543                                 // we have found one
1544                                 NewCellsList[NewNumberOfCellsInList]=CellsList[l1];
1545                                 NewNumberOfCellsInList++;
1546                                 break;
1547                               }
1548                             }
1549                         NumberOfCellsInList = NewNumberOfCellsInList;
1550
1551                         delete [] CellsList;
1552                         CellsList = NewCellsList;
1553                       }
1554
1555                       if (NumberOfCellsInList > 0) { // We have found some elements !
1556                         int CellNumber = CellsList[0];
1557                         //delete [] CellsList;
1558                         if (NumberOfCellsInList>1)
1559                           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one other Cell ("<<NumberOfCellsInList<<") !"));
1560
1561                         if (NumberOfCellsInList==1)
1562                           {
1563                             ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=CellNumber;
1564
1565                             // we search sub cell number in this cell to not calculate it another time
1566                             // which type ?
1567                             CELLMODEL Type2;
1568                             for (int l=0; l<_numberOfTypes; l++)
1569                               if (CellNumber < _count[l+1]) {
1570                                 Type2=_type[l];
1571                                 break;
1572                               }
1573                             //int sub_cell_count2 = Type2.get_entities_count(1);
1574                             //int nodes_cell_count2 = Type2.get_nodes_count();
1575                             bool find2 = false;
1576                             for (int l=1; l<=Type2.getNumberOfConstituents(1);l++) { // on all sub cell
1577                               int counter = 0;
1578                               for (int m=1; m<=Type2.getConstituentType(1,l)%100; m++)
1579                                 for (int n=1; n<=Type.getConstituentType(1,k)%100; n++)
1580                                   {
1581                                     if (_nodal->getIJ(CellNumber,Type2.getNodeConstituent(1,l,m)) == NodesLists[n-1])
1582                                       counter++;
1583                                   }
1584                               if (counter==Type.getConstituentType(1,k)%100)
1585                                 {
1586                                   descend_connectivity[descend_connectivity_index[CellNumber-1]+l-2]=-1*TotalNumberOfSubCell; // because, see it in other side !
1587                                   find2 = true;
1588                                 }
1589                               if (find2)
1590                                 break;
1591                             }
1592                             if (!find2)
1593                               MESSAGE_MED(LOC<<"ERROR ERROR ERROR ERROR ERROR : we find any subcell !!!"); // exception ?
1594                           }
1595                       } else {
1596                         ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=0;
1597                       }
1598                       delete [] CellsList;
1599                     }
1600
1601                   delete [] NodesLists;
1602                 }
1603               }
1604         }
1605       // we adjust _constituent
1606       int NumberOfConstituent=0;
1607       int SizeOfConstituentNodal=0;
1608       for (int i=0;i<_constituent->_numberOfTypes; i++) {
1609         NumberOfConstituent += tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1;
1610         SizeOfConstituentNodal += (tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1)*_constituent->_type[i].getNumberOfNodes();
1611       }
1612       // we built new _nodal attribute in _constituent
1613       //MEDSKYLINEARRAY *ConstituentNodal = new MEDSKYLINEARRAY(NumberOfConstituent,SizeOfConstituentNodal);
1614       //const int *ConstituentNodalValue = ConstituentNodal->getValue();
1615       //const int *ConstituentNodalIndex = ConstituentNodal->getIndex();
1616       int *ConstituentNodalValue = new int[SizeOfConstituentNodal];
1617       int *ConstituentNodalIndex = new int[NumberOfConstituent+1];
1618       ConstituentNodalIndex[0]=1;
1619       // we build _reverseDescendingConnectivity
1620       //_reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,2*NumberOfConstituent);
1621       //const int *reverseDescendingConnectivityValue = _reverseDescendingConnectivity->getValue();
1622       //const int *reverseDescendingConnectivityIndex = _reverseDescendingConnectivity->getIndex();
1623       int *reverseDescendingConnectivityValue = new int[2*NumberOfConstituent];
1624       int *reverseDescendingConnectivityIndex = new int[NumberOfConstituent+1];
1625       reverseDescendingConnectivityIndex[0]=1;
1626
1627       // first constituent type
1628       for(int j=0; j<tmp_NumberOfConstituentsForeachType[0]; j++) {
1629         ConstituentNodalIndex[j+1]=ConstituentNodalConnectivityIndex[j+1];
1630         for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1631           ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[k-1];
1632         }
1633         reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1634         for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1635           reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[k-1];
1636         }
1637       }
1638       // second type if any
1639       if (_constituent->_numberOfTypes==2) {
1640         int offset = _constituent->_count[1]-tmp_NumberOfConstituentsForeachType[0]-1;
1641         int offset1=offset*_constituent->_type[0].getNumberOfNodes();
1642         int offset2=offset*2;
1643         int NumberOfNodesPerConstituent = _constituent->_type[1].getNumberOfNodes();
1644         for(int j=tmp_NumberOfConstituentsForeachType[0]; j<NumberOfConstituent; j++) {
1645           ConstituentNodalIndex[j+1]=ConstituentNodalIndex[j]+NumberOfNodesPerConstituent;
1646           for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1647             ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[offset1+k-1];
1648           }
1649           reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1650           for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1651             reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[offset2+k-1];
1652           }
1653         }
1654         _constituent->_count[2]=NumberOfConstituent+1;
1655         // we correct _descending to adjust face number
1656         for(int j=0;j<DescendingSize;j++)
1657           if (abs(descend_connectivity[j])>tmp_NumberOfConstituentsForeachType[0]) {
1658             if ( descend_connectivity[j] > 0 )
1659               descend_connectivity[j]-=offset;
1660             else
1661               descend_connectivity[j]+=offset;
1662           }
1663       }
1664
1665       delete [] ConstituentNodalConnectivityIndex;
1666       delete [] ConstituentNodalConnectivity;
1667       delete [] ReverseDescendingConnectivityValue;
1668       if (_constituent->_numberOfTypes > 0)
1669         _constituent->_count[1]=tmp_NumberOfConstituentsForeachType[0]+1;
1670       delete [] tmp_NumberOfConstituentsForeachType;
1671
1672       _descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,
1673                                         DescendingSize,
1674                                         descend_connectivity_index,
1675                                         descend_connectivity);
1676       delete [] descend_connectivity_index;
1677       delete [] descend_connectivity;
1678
1679       ////
1680       vector<int> PolyDescending;
1681       vector<int> Reversedescendingconnectivityvalue(reverseDescendingConnectivityValue,reverseDescendingConnectivityValue + 2*NumberOfConstituent);
1682       vector<int> Reversedescendingconnectivityindex(reverseDescendingConnectivityIndex,reverseDescendingConnectivityIndex + NumberOfConstituent);
1683       delete [] reverseDescendingConnectivityValue;
1684       delete [] reverseDescendingConnectivityIndex;
1685
1686
1687       // polygons (2D mesh)
1688
1689       vector<int> Constituentnodalvalue(ConstituentNodalValue,ConstituentNodalValue + SizeOfConstituentNodal);
1690       vector<int> Constituentnodalindex(ConstituentNodalIndex,ConstituentNodalIndex + NumberOfConstituent+1);
1691       delete [] ConstituentNodalValue;
1692       delete [] ConstituentNodalIndex;
1693       int NumberOfNewSeg = 0;
1694
1695       for (int i=0; i <getNumberOfPolygons(); i++) // for each polygon
1696         {
1697           const int * vector_begin = &_polygonsNodal->getValue()[_polygonsNodal->getIndex()[i]-1];
1698           int vector_size = _polygonsNodal->getIndex()[i+1]-_polygonsNodal->getIndex()[i]+1;
1699           vector<int> myPolygon(vector_begin, vector_begin+vector_size);
1700           myPolygon[myPolygon.size()-1] = myPolygon[0]; // because first and last point make a segment
1701
1702           for (int j=0; j<(int)myPolygon.size()-1; j++) // for each segment of polygon
1703             {
1704               MEDMODULUSARRAY segment_poly(2,&myPolygon[j]);
1705               int ret_compare = 0;
1706
1707               // we search it in existing segments
1708
1709               for (int k=0; k<(int)Constituentnodalindex.size()-1; k++)
1710                 {
1711                   MEDMODULUSARRAY segment(2,&Constituentnodalvalue[0] + Constituentnodalindex[k]-1);
1712                   ret_compare = segment_poly.compare(segment);
1713                   if (ret_compare) // segment_poly already exists
1714                     {
1715                       PolyDescending.push_back(ret_compare*(k+1)); // we had it to the connectivity
1716                       insert_vector(Reversedescendingconnectivityvalue, 2*k+1, i+1 + getNumberOf(MED_CELL,MED_ALL_ELEMENTS)); // add polygon i to reverse descending connectivity for segment_poly (in 2sd place)
1717                       break;
1718                     }
1719                 }
1720
1721               // segment_poly must be created
1722
1723               if (!ret_compare)
1724                 {
1725                   NumberOfNewSeg++;
1726                   PolyDescending.push_back(NumberOfConstituent+NumberOfNewSeg); // we had it to the connectivity
1727                   Constituentnodalvalue.push_back(segment_poly[0]);
1728                   Constituentnodalvalue.push_back(segment_poly[1]);
1729                   Constituentnodalindex.push_back(Constituentnodalindex[Constituentnodalindex.size()-1] + 2); // we have only segments
1730                   insert_vector(Reversedescendingconnectivityvalue, 2*(NumberOfConstituent+NumberOfNewSeg-1), i+1 + getNumberOf(MED_CELL,MED_ALL_ELEMENTS)); // add polygon i to reverse descending connectivity for segment_poly (in 1st place)
1731                   insert_vector(Reversedescendingconnectivityindex, NumberOfConstituent+NumberOfNewSeg-1, 2*(NumberOfConstituent+NumberOfNewSeg-1)+1); // idem with index
1732                 }
1733             }
1734         }
1735
1736       if (getNumberOfPolygons() > 0)
1737         {
1738           _polygonsDescending = new MEDSKYLINEARRAY(getNumberOfPolygons(),_polygonsNodal->getLength(),_polygonsNodal->getIndex(),&PolyDescending[0]); // index are the same for polygons nodal and descending connectivities
1739
1740           NumberOfConstituent += NumberOfNewSeg;
1741           SizeOfConstituentNodal += 2*NumberOfNewSeg;
1742           _constituent->_count[1] = NumberOfConstituent+1;
1743         }
1744
1745
1746       // polyhedron (3D mesh)
1747
1748       vector<int> Constituentpolygonsnodalvalue;
1749       vector<int> Constituentpolygonsnodalindex(1,1);
1750       int NumberOfNewFaces = 0; // by convention new faces are polygons
1751       //offset to switch between all types and classical types.
1752       int offsetCell = getNumberOf(MED_CELL, MED_ALL_ELEMENTS);
1753       int *tabRes = new int[1000]; //temporay array for intersection calculation
1754
1755       for (int i=0; i<getNumberOfPolyhedron(); i++) // for each polyhedron
1756         {
1757           // we create a POLYHEDRONARRAY containing only this polyhedra
1758           int myNumberOfFaces = _polyhedronNodal->getPolyhedronIndex()[i+1]-_polyhedronNodal->getPolyhedronIndex()[i];
1759           int myNumberOfNodes = _polyhedronNodal->getFacesIndex()[_polyhedronNodal->getPolyhedronIndex()[i+1]-1]-_polyhedronNodal->getFacesIndex()[_polyhedronNodal->getPolyhedronIndex()[i]-1];
1760           POLYHEDRONARRAY myPolyhedra(1,myNumberOfFaces,myNumberOfNodes);
1761 //CCRT    vector<int> myFacesIndex(_polyhedronNodal->getFacesIndex() + _polyhedronNodal->getPolyhedronIndex()[i]-1, _polyhedronNodal->getFacesIndex() + _polyhedronNodal->getPolyhedronIndex()[i]-1 + myNumberOfFaces+1);
1762           vector<med_int> myFacesIndex(_polyhedronNodal->getFacesIndex() + _polyhedronNodal->getPolyhedronIndex()[i]-1, _polyhedronNodal->getFacesIndex() + _polyhedronNodal->getPolyhedronIndex()[i]-1 + myNumberOfFaces+1);
1763           for (int j=myNumberOfFaces; j>=0; j--)
1764             myFacesIndex[j] -= myFacesIndex[0]-1;
1765           myPolyhedra.setFacesIndex(&myFacesIndex[0]);
1766           myPolyhedra.setNodes(_polyhedronNodal->getNodes() + _polyhedronNodal->getFacesIndex()[_polyhedronNodal->getPolyhedronIndex()[i]-1]-1);
1767
1768           for (int j=0; j<myPolyhedra.getNumberOfFaces(); j++) // for each face of polyhedra
1769             {
1770               int myFaceNumberOfNodes = myPolyhedra.getFacesIndex()[j+1]-myPolyhedra.getFacesIndex()[j];
1771 //CCRT
1772 #if defined(IRIX64) || defined(OSF1) || defined(VPP5000) || defined(PCLINUX64)
1773               int ii ;
1774               const med_int * Nodes = myPolyhedra.getNodes() ;
1775               int * tmp_Nodes = new int[myPolyhedra.getNumberOfNodes()] ;
1776               for ( ii = 0 ; ii < myPolyhedra.getNumberOfNodes() ; ii++ )
1777                  tmp_Nodes[ii] = Nodes[ii] ;
1778               const med_int * FacesIndex = myPolyhedra.getFacesIndex() ;
1779               int * tmp_FacesIndex = new int[myPolyhedra.getNumberOfFaces()+1] ;
1780               for ( ii = 0 ; ii < myPolyhedra.getNumberOfFaces() ; ii++ )
1781                  tmp_FacesIndex[ii] = FacesIndex[ii] ;
1782 //CCRT : copy of Nodes
1783               MEDMODULUSARRAY face_poly(myFaceNumberOfNodes,tmp_Nodes + tmp_FacesIndex[j]-1);
1784 //CCRT              delete [] tmp_Nodes ;
1785               delete [] tmp_FacesIndex ;
1786 #else
1787               MEDMODULUSARRAY face_poly(myFaceNumberOfNodes,myPolyhedra.getNodes() + myPolyhedra.getFacesIndex()[j]-1);
1788 #endif
1789               int ret_compare = 0;
1790
1791               // we search it in existing faces
1792
1793               // we search first in TRIA3 and QUAD4
1794               if (myFaceNumberOfNodes == 3 || myFaceNumberOfNodes == 4)
1795                 {
1796                   int Begin = -1; // first TRIA3 or QUAD4
1797                   int Number = 0; // number of TRIA3 or QUAD4
1798                   for (int k=0; k<(int)Constituentnodalindex.size()-1; k++)
1799                     if (Constituentnodalindex[k+1]-Constituentnodalindex[k] == myFaceNumberOfNodes)
1800                       {
1801                         if (Begin == -1)
1802                           Begin = k;
1803                         Number++;
1804                       }
1805
1806                   for (int k=0; k<Number; k++)
1807                     {
1808                       MEDMODULUSARRAY face(myFaceNumberOfNodes,&Constituentnodalvalue[0] + Constituentnodalindex[Begin+k]-1);
1809                       ret_compare = face_poly.compare(face);
1810                       if (ret_compare)
1811                         {
1812                           PolyDescending.push_back(ret_compare*(Begin+k+1)); // we had it to the connectivity
1813                           insert_vector(Reversedescendingconnectivityvalue, 2*(Begin+k)+1, i+1 + getNumberOf(MED_CELL,MED_ALL_ELEMENTS)); // add polyhedra i to reverse descending connectivity for face_poly (in 2sd place)
1814                           break;
1815                         }
1816                     }
1817                 }
1818
1819               // we search last in POLYGONS
1820               if (!ret_compare)
1821               {
1822                 int lgth;
1823                 const int *facePolyTab=face_poly.getArray(lgth);
1824                 int nbOfCandidatesCell = ReverseNodalConnectivityIndex[facePolyTab[0]] -
1825                                          ReverseNodalConnectivityIndex[facePolyTab[0]-1];
1826                 const int *candidatesCell = ReverseNodalConnectivityValue +
1827                                             ReverseNodalConnectivityIndex[facePolyTab[0]-1] - 1;
1828                 memcpy(tabRes,candidatesCell,nbOfCandidatesCell*sizeof(int));
1829                 int lgth2=nbOfCandidatesCell;
1830                 for (int k=1;k<lgth && lgth2!=0;k++)
1831                 {
1832                   nbOfCandidatesCell = ReverseNodalConnectivityIndex[facePolyTab[k]] -
1833                                        ReverseNodalConnectivityIndex[facePolyTab[k]-1];
1834                   candidatesCell = ReverseNodalConnectivityValue +
1835                                    ReverseNodalConnectivityIndex[facePolyTab[k]-1] - 1;
1836                   mergeOrderedTabs(tabRes,lgth2,candidatesCell,nbOfCandidatesCell,tabRes,lgth2);
1837                 }
1838                 if (lgth2<=1)
1839                   ret_compare=0;//here normally tabRes[0]==offsetCell+i+1
1840                 else //> 2 should never happend : A face is shared by more than 2 polyhedrons...
1841                 {
1842                   if (tabRes[0] == offsetCell+i+1) //as tabRes is ordered by construction tabRes[1] > tabRes[0] so the current
1843                     // face is shared with an another cell whose id > current id. So let's create
1844                     ret_compare=0;
1845                   else
1846                   {//tabRes[0]<Constituentpolygonsnodalindex.size()-1 that is to say the current face has been built previously.
1847                     const int *facesConstitutingAlreadyBuiltPolyh = &PolyDescending[0] + _polyhedronNodal->getPolyhedronIndex()[tabRes[0]-offsetCell-1] - 1;
1848                     int nbOfFacesConstitutingAlreadyBuiltPolyh = _polyhedronNodal->getPolyhedronIndex()[tabRes[0]-offsetCell] - _polyhedronNodal->getPolyhedronIndex()[tabRes[0]-offsetCell-1];
1849                     for (int k1=0; k1<nbOfFacesConstitutingAlreadyBuiltPolyh && (ret_compare==0); k1++)
1850                     {
1851                       int curFaceId=facesConstitutingAlreadyBuiltPolyh[k1];
1852                       if(curFaceId>NumberOfConstituent)//In other case it is not a polyhedron : no chance to fit if you see comment 30 lines behind.
1853                       {
1854                         int nbOfNodesForCurrentFace =
1855                           Constituentpolygonsnodalindex[curFaceId-NumberOfConstituent]
1856                           - Constituentpolygonsnodalindex[curFaceId-NumberOfConstituent-1];
1857                         MEDMODULUSARRAY face (nbOfNodesForCurrentFace,&Constituentpolygonsnodalvalue[0]+
1858                                               Constituentpolygonsnodalindex[curFaceId-NumberOfConstituent-1]-1);
1859                         ret_compare = face_poly.compare(face);
1860                         if (ret_compare)
1861                         {
1862                           PolyDescending.push_back(ret_compare*curFaceId); // we had it to the connectivity
1863                           insert_vector(Reversedescendingconnectivityvalue, 2*(curFaceId-1)+1,
1864                                         i + 1 + getNumberOf(MED_CELL,MED_ALL_ELEMENTS));
1865                         }
1866                       }
1867                     }
1868                   }
1869                 }
1870               }
1871
1872               // if not found, face_poly must be created
1873
1874               if (!ret_compare)
1875                 {
1876                   NumberOfNewFaces++;
1877                   PolyDescending.push_back(NumberOfConstituent+NumberOfNewFaces); // we had it to the connectivity
1878                   for (int k=0; k<myFaceNumberOfNodes; k++)
1879                     Constituentpolygonsnodalvalue.push_back(face_poly[k]);
1880                   Constituentpolygonsnodalindex.push_back(Constituentpolygonsnodalindex[Constituentpolygonsnodalindex.size()-1]+myFaceNumberOfNodes);
1881                   insert_vector(Reversedescendingconnectivityvalue, 2*(NumberOfConstituent+NumberOfNewFaces-1), i+1 + getNumberOf(MED_CELL,MED_ALL_ELEMENTS)); // add polyhedra i to reverse descending connectivity for face_poly (in 1st place)
1882                   insert_vector(Reversedescendingconnectivityindex, NumberOfConstituent+NumberOfNewFaces-1, 2*(NumberOfConstituent+NumberOfNewFaces-1)+1); // idem with index
1883                 }
1884             }
1885         }
1886       delete [] tabRes;
1887
1888       if (getNumberOfPolyhedron() > 0)
1889         {
1890 #if defined(IRIX64) || defined(OSF1) || defined(VPP5000) || defined(PCLINUX64)
1891           int * tmp_PolyhedronIndex = new int[getNumberOfPolyhedron()+1] ;
1892           const MED_EN::med_int * PolyhedronIndex = _polyhedronNodal->getPolyhedronIndex() ;
1893           int ii ;
1894           for ( ii = 0 ; ii < getNumberOfPolyhedron()+1 ; ii++ )
1895              tmp_PolyhedronIndex[ii] = PolyhedronIndex[ii] ;
1896 //CCRT : copy of PolyhedronIndex
1897           _polyhedronDescending = new MEDSKYLINEARRAY(getNumberOfPolyhedron(),_polyhedronNodal->getNumberOfFaces(),tmp_PolyhedronIndex,&PolyDescending[0]); // polyhedron index are the same for nodal and descending connectivities
1898 #else
1899           _polyhedronDescending = new MEDSKYLINEARRAY(getNumberOfPolyhedron(),_polyhedronNodal->getNumberOfFaces(),_polyhedronNodal->getPolyhedronIndex(),&PolyDescending[0]); // polyhedron index are the same for nodal and descending connectivities
1900 #endif
1901
1902           if (_constituent->_polygonsNodal != NULL)
1903             delete [] _constituent->_polygonsNodal;
1904           _constituent->_polygonsNodal = new MEDSKYLINEARRAY(Constituentpolygonsnodalindex.size()-1,Constituentpolygonsnodalvalue.size(),&Constituentpolygonsnodalindex[0],&Constituentpolygonsnodalvalue[0]);
1905         }
1906
1907       // delete _constituent->_nodal
1908       _constituent->_nodal = new MEDSKYLINEARRAY(NumberOfConstituent,
1909                                                  SizeOfConstituentNodal,
1910                                                  &Constituentnodalindex[0],
1911                                                  &Constituentnodalvalue[0]);
1912
1913       int NumberOfConstituentWithPolygons = NumberOfConstituent + NumberOfNewFaces;
1914       Reversedescendingconnectivityindex.push_back(Reversedescendingconnectivityindex[Reversedescendingconnectivityindex.size()-1]+2); // we complete the index
1915       _reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituentWithPolygons+1,
1916                                                            2*NumberOfConstituentWithPolygons,
1917                                                            &Reversedescendingconnectivityindex[0],
1918                                                            &Reversedescendingconnectivityvalue[0]);
1919       ////
1920       _isDescendingConnectivityPartial = false;
1921     }
1922   END_OF_MED(LOC);
1923   }
1924 }
1925
1926 void CONNECTIVITY::addToDescendingConnectivity(const set<int>& nodes,
1927                                                multimap<int,int>& descending,
1928                                                int iglobal_cell,
1929                                                const CONNECTIVITY_HashMap & face_map )
1930 {
1931         int dimension = getEntityDimension();
1932         vector<int> signature (dimension);
1933         set<int>::const_iterator iter=nodes.begin();
1934         for (int i=0; i< dimension;i++)
1935                 {
1936                         signature[i]=*iter;
1937                         iter++;
1938                 }
1939
1940
1941         CONNECTIVITY_HashMap::const_iterator itermap=face_map.find(signature);
1942         CONNECTIVITY_HashMap::const_iterator iterend=face_map.end();
1943
1944
1945         if (itermap!=iterend)
1946                 descending.insert(make_pair(iglobal_cell,itermap->second));
1947
1948 }
1949
1950 /*! This method calculates the descending connectivity without creating missing elements. It only maps the constituent elements that are described in the nodal representation.
1951 For instance, let us consider the following mesh with no MED_EDGE elements on the inner edges.
1952
1953 +----1----+----2----+
1954 |         |         |
1955 8    1    |    2    3
1956 |         |         |
1957 +---------+---------+
1958 |         |         |
1959 7    3    |    4    4
1960 |         |         |
1961 +----6----+----5----+
1962
1963
1964 calculatePartialDescendingConnectivity()
1965 will return
1966 1:1,8
1967 2:2,3
1968 3:6,7
1969 4:4,5
1970
1971 whereas calculateDescendingConnectivity()
1972 will create new edges, renumbering existing ones and return
1973 1:1,2,3,4
1974 2:5,6,7,2
1975 3:3,8,9,10
1976 4:7,11,12,8
1977
1978 +----1----+----5----+
1979 |         |         |
1980 4    1    2    2    6
1981 |         |         |
1982 +----3----+----7----+
1983 |         |         |
1984 10    3   8    4    11
1985 |         |         |
1986 +----9----+----12---+
1987 */
1988
1989 void CONNECTIVITY::calculatePartialDescendingConnectivity()
1990 {
1991         ////////////////////////////////////////////////////////////////////////////
1992         // First stage : creating hash_map referencing faces with the triplet
1993         // of the lowest order nodes as a key and the global face number as a value
1994
1995         CONNECTIVITY_HashMap face_map;
1996         const medGeometryElement* face_types = _constituent->getGeometricTypes(_constituent->getEntity());
1997         int nbface_types=_constituent->getNumberOfTypes(_constituent->getEntity());
1998         int iglobal_face=1;
1999
2000         int dimension = getEntityDimension();
2001
2002         //classical types
2003         for (int itype=0; itype<nbface_types; itype++)
2004                 {
2005                         medGeometryElement face_type=face_types[itype];
2006                         int nbnodes=face_type%100;
2007                         int nbfaces=_constituent->getNumberOf(_constituent->getEntity(),face_type);
2008                         for (int iface=0;iface<nbfaces;iface++)
2009                                 {
2010                                         set<int> nodes;
2011                                         const int* index=_constituent->_nodal->getIndex();
2012                                         const int* conn=_constituent->_nodal->getValue();
2013                                         for (int inode=0;inode<nbnodes;inode++)
2014                                                 nodes.insert(conn[index[iglobal_face-1]-1+inode]);
2015                                         vector<int> signature (dimension);
2016                                         set<int>::iterator iter=nodes.begin();
2017                                         for (int i=0; i< dimension;i++)
2018                                                 {
2019                                                         signature[i]=*iter;
2020                                                         iter++;
2021                                                 }
2022                                         face_map.insert(make_pair(signature,iglobal_face));
2023                                         iglobal_face++;
2024                                 }
2025                 }
2026
2027         //polygons
2028         if (dimension==3 && _polygonsNodal!=0)
2029                 {
2030                         int nbfaces=_constituent->getNumberOfPolygons();
2031                         for (int iface=0;iface<nbfaces;iface++)
2032                                 {
2033                                         set<int> nodes;
2034                                         const int* index=_constituent->_polygonsNodal->getIndex();
2035                                         int nbnodes=index[iface+1]-index[iface];
2036                                         const int* conn=_constituent->_polygonsNodal->getValue();
2037                                         for (int inode=0;inode<nbnodes;inode++)
2038                                                 nodes.insert(conn[index[iface]-1+inode]);
2039                                         vector<int> signature (dimension);
2040                                         set<int>::iterator iter=nodes.begin();
2041                                         for (int i=0; i< dimension;i++)
2042                                                 {
2043                                                         signature[i]=*iter;
2044                                                         iter++;
2045                                                 }
2046                                         face_map.insert(make_pair(signature,iglobal_face));
2047                                         iglobal_face++;
2048                                 }
2049                 }
2050
2051         ////////////////////////////////////////////////////////////////////////////
2052         //Second stage : going through all the faces of the cells and
2053         // connecting them to the hash_map created in the first stage
2054
2055         multimap<int,int> descending; //map storing the descending connectivity
2056
2057         //classical types
2058         int nbcell_types=getNumberOfTypes(_entity);
2059         const medGeometryElement* cell_types=getGeometricTypes(_entity);
2060         int iglobal_cell=1;
2061         for (int itype=0; itype<nbcell_types; itype++)
2062                 {
2063                         medGeometryElement cell_type=cell_types[itype];
2064                         int nbcells=getNumberOf(_entity,cell_type);
2065
2066                         CELLMODEL cellmodel=CELLMODEL_Map::retrieveCellModel(cell_type);
2067                         for (int icell=0;icell<nbcells;icell++)
2068                                 {
2069                                         int nbfaces=cellmodel.getNumberOfConstituents(1);
2070                                         for (int iface=0; iface<nbfaces;iface++)
2071                                                 {
2072                                                         set<int> nodes;
2073                                                         const int* local_index=cellmodel.getNodesConstituent(1,iface+1);
2074                                                         medGeometryElement face_type = cellmodel.getConstituentType(1,iface+1);
2075                                                         const int* index=_nodal->getIndex();
2076                                                         const int* conn=_nodal->getValue();
2077                                                         int nbnodes=face_type%100;
2078                                                         for (int inode=0;inode<nbnodes;inode++)
2079                                                                 nodes.insert(conn[index[iglobal_cell-1]-1+local_index[inode]-1]);
2080                                                         addToDescendingConnectivity(nodes,descending,iglobal_cell,face_map);
2081                                                 }
2082                                         iglobal_cell++;
2083                                 }
2084                 }
2085
2086         // polygon case
2087         if (dimension==2 && _polygonsNodal!=0)
2088         {
2089           int nbcells=getNumberOfElementsWithPoly(_entity,MED_POLYGON);
2090           const int* index = _polygonsNodal->getIndex();
2091           const int* conn = _polygonsNodal->getValue();
2092
2093           for (int icell=0;icell<nbcells;icell++)
2094           {
2095             int nbnodes=index[icell+1]-index[icell];
2096             int nbfaces=nbnodes;
2097             for (int iface=0; iface<nbfaces;iface++)
2098             {
2099               set<int> nodes;
2100               if (iface+1!=nbfaces)
2101               {
2102                 nodes.insert(conn[index[icell]-1+iface]);
2103                 nodes.insert(conn[index[icell]-1+iface+1]);
2104               }
2105               else
2106               {
2107                 nodes.insert(conn[index[icell]-1+iface]);
2108                 nodes.insert(conn[index[icell]-1]);
2109               }
2110               addToDescendingConnectivity(nodes,descending,iglobal_cell,face_map);
2111             }
2112             iglobal_cell++;
2113           }
2114         }
2115
2116         //polyhedra case
2117         if (dimension==3 && _polyhedronNodal!=0)
2118         {
2119           int nbcells = getNumberOfElementsWithPoly(_entity,MED_POLYHEDRA);
2120           const MED_EN::med_int* index = _polyhedronNodal->getPolyhedronIndex();
2121           const MED_EN::med_int* conn = _polyhedronNodal->getNodes();
2122           const MED_EN::med_int* face_index = _polyhedronNodal->getFacesIndex();
2123
2124           for (int icell = 0; icell < nbcells; icell++)
2125           {
2126             for (int iface = face_index[icell]; iface < face_index[icell+1]; iface++)
2127             {
2128               set<int> nodes;
2129               for (int inode = index[iface-1]; iface < index[iface]; iface++)
2130                 nodes.insert(conn[inode-1]);
2131               addToDescendingConnectivity(nodes,descending,iglobal_cell,face_map);
2132             }
2133             iglobal_cell++;
2134           }
2135         }
2136
2137
2138         ////////////////////////////////////////////////////////////////////////////
2139         // Third stage : reorganizing the descending data to store it in a medskylinearray
2140
2141         vector<int> index;
2142         vector<int> value;
2143         index.push_back(1);
2144
2145         //the number of cells is given by the number
2146         //obtained by browsing all the cell types
2147         int nb_cells = iglobal_cell-1;
2148
2149         for (int icell=0; icell< nb_cells;icell++)
2150                 {
2151                         multimap<int,int>::iterator beginning_of_range = descending.lower_bound(icell);
2152                         multimap<int,int>::iterator end_of_range = descending.upper_bound(icell);
2153                         int nb=0;
2154                         for (multimap<int,int>::iterator iter2 = beginning_of_range; iter2!=end_of_range;iter2++)
2155                                 {
2156                                         value.push_back(iter2->second);
2157                                         nb++;
2158                                 }
2159                         index.push_back(index[index.size()-1]+nb);
2160                 }
2161         _descending=new MEDSKYLINEARRAY(index.size()-1,value.size(),&index[0],&value[0]);
2162         _isDescendingConnectivityPartial = true;
2163 }
2164
2165
2166 /*! Not implemented yet */
2167 //--------------------------------------------------------------------//
2168 void CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity)
2169 //--------------------------------------------------------------------//
2170 {
2171   const char* LOC = "CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity) : ";
2172   BEGIN_OF_MED(LOC);
2173
2174   MESSAGE_MED(PREFIX_MED<<"method not yet implemented " << myConnectivity._entity);
2175   // Mesh dimension !
2176
2177   END_OF_MED(LOC);
2178   return;
2179 }
2180
2181
2182 /*!
2183   Returns the geometry of an element given by its entity type & its global number.
2184
2185   Example : medGeometryElement myType = myConnectivity.getElementType(MED_CELL,35);
2186 */
2187 //--------------------------------------------------------------------//
2188 medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity,int globalNumber) const
2189 //--------------------------------------------------------------------//
2190 {
2191   const char * LOC = "medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity, int globalNumber) const : ";
2192   BEGIN_OF_MED(LOC);
2193   int globalNumberMin = 1;
2194   int globalNumberMax ;
2195
2196   if (_entity==Entity) globalNumberMax = _count[_numberOfTypes];
2197   else if (_constituent!=NULL) globalNumberMax = _constituent->_count[_constituent->_numberOfTypes];
2198   else
2199     throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
2200
2201   // The globalNumber must verify  : 1 <=  globalNumber < _count[_numberOfTypes] (== totalNumberOfElement+1)
2202
2203   if ( (globalNumber < globalNumberMin) || (globalNumber >  globalNumberMax-1 )  )
2204     throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "globalNumber must be between >= |"
2205                                  << globalNumberMin <<"| and <= |" << globalNumberMax-1 ));
2206
2207   if (_entity==Entity) {
2208     for(int i=1; i<=_numberOfTypes;i++)
2209       if (globalNumber<_count[i])
2210         return _geometricTypes[i-1];
2211   }
2212   else if (_constituent!=NULL)
2213     return _constituent->getElementType(Entity,globalNumber);
2214   else
2215     throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
2216   throw MEDEXCEPTION(LOCALIZED("getElementType : Wrong Number !"));
2217
2218   END_OF_MED(LOC);
2219 }
2220
2221 /*!
2222   Idem as getElementType method except that it manages polygon and polyhedron.
2223 */
2224 //--------------------------------------------------------------------//
2225 medGeometryElement CONNECTIVITY::getElementTypeWithPoly(medEntityMesh Entity,int globalNumber) const
2226 //--------------------------------------------------------------------//
2227 {
2228   int globalNumberMaxOfClassicType;
2229   if(_entity==Entity)
2230     {
2231       if (_count==0)
2232         globalNumberMaxOfClassicType = 1;
2233       else
2234         globalNumberMaxOfClassicType=_count[_numberOfTypes];
2235       if(globalNumber>=1)
2236         {
2237           if(globalNumber<globalNumberMaxOfClassicType)
2238             {
2239               for(int i=1; i<=_numberOfTypes;i++)
2240                 if (globalNumber<_count[i])
2241                   return _geometricTypes[i-1];
2242               throw MEDEXCEPTION("never happens just for compilo");
2243             }
2244           else
2245             {
2246               int localNumberInPolyArray=globalNumber-globalNumberMaxOfClassicType;
2247               int nbOfPol=getNumberOfElementOfPolyType(_entity);
2248               if(localNumberInPolyArray<nbOfPol)
2249                 return getPolyTypeRelativeTo();
2250               else
2251                 throw MEDEXCEPTION("getElementTypeWithPoly : unexisting element type");
2252             }
2253         }
2254       else
2255         throw MEDEXCEPTION("getElementTypeWithPoly : globalNumber < 1");
2256     }
2257   else
2258     {
2259       if(_constituent!=NULL)
2260         return _constituent->getElementTypeWithPoly(Entity,globalNumber);
2261       else
2262         throw MEDEXCEPTION("getElementTypeWithPoly : Entity not defined !");
2263     }
2264 }
2265
2266 ostream & MEDMEM::operator<<(ostream &os, CONNECTIVITY &co)
2267 {
2268     os << endl << "------------- Entity = ";
2269     switch (co._entity)
2270     {
2271         case MED_CELL:
2272             os << "MED_CELL";
2273             break;
2274         case MED_FACE:
2275             os << "MED_FACE";
2276             break;
2277         case MED_EDGE:
2278             os << "MED_EDGE";
2279             break;
2280         case MED_NODE:
2281             os << "MED_NODE";
2282             break;
2283         case MED_ALL_ENTITIES:
2284             os << "MED_ALL_ENTITIES";
2285             break;
2286         default:
2287             os << "Unknown";
2288             break;
2289     }
2290     os  << " -------------\n\nMedConnectivity : ";
2291     switch (co._typeConnectivity)
2292     {
2293         case MED_NODAL:
2294             os << "MED_NODAL\n";
2295             break;
2296         case MED_DESCENDING:
2297             os << "MED_DESCENDING\n";
2298             break;
2299         default:
2300             break;
2301     }
2302     os << "Entity dimension : " << co._entityDimension << endl;
2303     os << "Number of nodes : " << co._numberOfNodes << endl;
2304     os << "Number of types : " << co._numberOfTypes << endl;
2305     for (int i=0; i!=co._numberOfTypes ; ++i)
2306         os << "  -> Type " << co._geometricTypes[i] << " (" << co._type[i].getName() << ") : "
2307            << co._count[i+1]-co._count[i] << " elements" << endl;
2308
2309     if (co._typeConnectivity == MED_NODAL )
2310     {
2311         for (int i=0; i<co._numberOfTypes; i++)
2312         {
2313             os << endl << co._type[i].getName() << " : " << endl;
2314             int numberofelements = co._count[i+1]-co._count[i];
2315             const int * connectivity = co.getConnectivity(co._typeConnectivity, co._entity, co._geometricTypes[i]);
2316             int numberofnodespercell = co._geometricTypes[i]%100;
2317             for (int j=0;j<numberofelements;j++)
2318             {
2319                 os << setw(6) << j+1 << " : ";
2320                 for (int k=0;k<numberofnodespercell;k++)
2321                     os << connectivity[j*numberofnodespercell+k]<<" ";
2322                 os << endl;
2323             }
2324         }
2325
2326         os << endl << "MED_POLYGON : " << co.getNumberOfPolygons() << " polygons" << endl;
2327         if (co.getNumberOfPolygons() > 0)
2328           {
2329             const int* polygonsconnectivity = co.getPolygonsConnectivity(MED_NODAL,co.getEntity());
2330             const int* polygonsconnectivityindex = co.getPolygonsConnectivityIndex(MED_NODAL,co.getEntity());
2331
2332             for (int i=0; i<co.getNumberOfPolygons(); i++)
2333               {
2334                 int numberofnodesperpolygon = polygonsconnectivityindex[i+1]-polygonsconnectivityindex[i];
2335                 for (int j=0; j<numberofnodesperpolygon; j++)
2336                   os << polygonsconnectivity[polygonsconnectivityindex[i]-1+j] << " ";
2337                 os << endl;
2338               }
2339           }
2340
2341         os << endl << "MED_POLYHEDRA : " << co.getNumberOfPolyhedron() << " polyhedron" << endl;
2342         if (co.getNumberOfPolyhedron() > 0)
2343           {
2344             const int* polyhedronconnectivity = co.getPolyhedronConnectivity(MED_NODAL);
2345             const int* polyhedronfacesindex = co.getPolyhedronFacesIndex();
2346             const int* polyhedronindex = co.getPolyhedronIndex(MED_NODAL);
2347
2348             for (int i=0; i<co.getNumberOfPolyhedron(); i++)
2349               {
2350                 int numberoffacesperpolyhedra = polyhedronindex[i+1]-polyhedronindex[i];
2351                 for (int j=0; j<numberoffacesperpolyhedra; j++)
2352                   {
2353                     int numberofnodesperface = polyhedronfacesindex[polyhedronindex[i]-1+j+1]-polyhedronfacesindex[polyhedronindex[i]-1+j];
2354                     for (int k=0; k<numberofnodesperface; k++)
2355                       os << polyhedronconnectivity[polyhedronfacesindex[polyhedronindex[i]-1+j]-1+k] << " ";
2356                     if (j != numberoffacesperpolyhedra-1) os << "| ";
2357                   }
2358                 os << endl;
2359               }
2360           }
2361
2362     }
2363     else if (co._typeConnectivity == MED_DESCENDING)
2364     {
2365         int numberofelements = co.getNumberOf( co._entity , MED_ALL_ELEMENTS);
2366         if (numberofelements > 0)
2367           {
2368             const int *connectivity =  co.getConnectivity( co._typeConnectivity, co._entity, MED_ALL_ELEMENTS);
2369             const int *connectivity_index =  co.getConnectivityIndex( co._typeConnectivity, co._entity );
2370
2371             for ( int j=0; j!=numberofelements; j++ )
2372               {
2373                 os << "Element " << j+1 << " : ";
2374                 for ( int k=connectivity_index[j]; k!=connectivity_index[j+1]; k++ )
2375                   os << connectivity[k-1] << " ";
2376                 os << endl;
2377               }
2378           }
2379
2380         os << endl << "MED_POLYGON : " << co.getNumberOfPolygons() << " polygons" << endl;
2381         if (co.getNumberOfPolygons() > 0)
2382           {
2383             const int* polygonsconnectivity = co.getPolygonsConnectivity(MED_DESCENDING,co.getEntity());
2384             const int* polygonsconnectivityindex = co.getPolygonsConnectivityIndex(MED_DESCENDING,co.getEntity());
2385
2386             for (int i=0; i<co.getNumberOfPolygons(); i++)
2387               {
2388                 int numberofedgesperpolygon = polygonsconnectivityindex[i+1]-polygonsconnectivityindex[i];
2389                 for (int j=0; j<numberofedgesperpolygon; j++)
2390                   os << polygonsconnectivity[polygonsconnectivityindex[i]-1+j] << " ";
2391                 os << endl;
2392               }
2393           }
2394
2395         os << endl << "MED_POLYHEDRA : " << co.getNumberOfPolyhedron() << " polyhedron" << endl;
2396         if (co.getNumberOfPolyhedron() > 0)
2397           {
2398             const int* polyhedronconnectivity = co.getPolyhedronConnectivity(MED_DESCENDING);
2399             const int* polyhedronindex = co.getPolyhedronIndex(MED_DESCENDING);
2400
2401             for (int i=0; i<co.getNumberOfPolyhedron(); i++)
2402               {
2403                 int numberoffacesperpolyhedra = polyhedronindex[i+1]-polyhedronindex[i];
2404                 for (int j=0; j<numberoffacesperpolyhedra; j++)
2405                   os << polyhedronconnectivity[polyhedronindex[i]-1+j] << " ";
2406                 os << endl;
2407               }
2408           }
2409
2410     }
2411
2412     if (co._constituent)
2413         os << endl << *co._constituent << endl;
2414
2415     return os;
2416 }
2417
2418 /*!
2419   method that adds to vector 'nodes' all the nodes of polyhedron with id 'polyhedronId'.
2420   WARNING the returned pointer should be deallocated. Returned nodes and polyhedronId are in form [1,...]
2421  */
2422 int *CONNECTIVITY::getNodesOfPolyhedron(int polyhedronId, int& lgthOfTab) const
2423 {
2424   const int *nodes=getPolyhedronConnectivity(MED_EN::MED_NODAL);
2425   const int *faces=getPolyhedronFacesIndex();
2426   const int *glob=getPolyhedronIndex(MED_EN::MED_NODAL);
2427   int offsetWithClassicType=getNumberOf(_entity,MED_ALL_ELEMENTS);
2428   set<int> retInSet;
2429   if (polyhedronId<offsetWithClassicType || polyhedronId> getNumberOfElementsWithPoly (MED_CELL, MED_ALL_ELEMENTS))
2430     throw MEDEXCEPTION("Polyhedron ID does not match a polyhedron in the element range");
2431   int startFace=glob[polyhedronId-offsetWithClassicType-1]-1;
2432   int endFace=glob[polyhedronId-offsetWithClassicType]-1;
2433   int i;
2434   for(i=startFace;i!=endFace;i++)
2435     {
2436       for(int j=faces[i];j<faces[i+1];j++)
2437         retInSet.insert(nodes[j-1]);
2438     }
2439   lgthOfTab=retInSet.size();
2440   int *ret=new int[lgthOfTab];
2441   set<int>::iterator iter=retInSet.begin();
2442   i=0;
2443   for(iter=retInSet.begin();iter!=retInSet.end();iter++)
2444     ret[i++]=*iter;
2445   return ret;
2446 }
2447
2448 /*!
2449   Idem as MESH::getNodesOfPolyhedron except that returned nodes are sorted by face. 'nbOfNodesPerFaces' is an array of size 'nbOfFaces'.
2450   Returned int** has a size of 'nbOfNodesPerFaces' too, and for each element j in int** the size is  nbOfNodesPerFaces[j].
2451   Warning both returned 'nbOfNodesPerFaces' and returned value should be deallocated. Returned nodes and 'polyhedronId' are in form [1,...]
2452  */
2453 int **CONNECTIVITY::getNodesPerFaceOfPolyhedron(int polyhedronId, int& nbOfFaces, int* & nbOfNodesPerFaces) const
2454 {
2455   int i;
2456   const int *nodes=getPolyhedronConnectivity(MED_EN::MED_NODAL);
2457   const int *faces=getPolyhedronFacesIndex();
2458   const int *glob=getPolyhedronIndex(MED_EN::MED_NODAL);
2459   int offsetWithClassicType=getNumberOf(_entity,MED_ALL_ELEMENTS);
2460   if (polyhedronId<offsetWithClassicType || polyhedronId> getNumberOfElementsWithPoly (MED_CELL, MED_ALL_ELEMENTS))
2461     throw MEDEXCEPTION("Polyhedron ID does not match a polyhedron in the element range");
2462
2463   int startFace=glob[polyhedronId-offsetWithClassicType-1]-1;
2464   nbOfFaces=glob[polyhedronId-offsetWithClassicType]-startFace-1;
2465   nbOfNodesPerFaces=new int[nbOfFaces];
2466   int **ret=new int* [nbOfFaces];
2467   for(i=0;i<nbOfFaces;i++)
2468     {
2469       int startNode=faces[startFace+i]-1;
2470       nbOfNodesPerFaces[i]=faces[startFace+i+1]-startNode-1;
2471       ret[i]=(int *)(nodes)+startNode;
2472     }
2473   return ret;
2474 }
2475
2476 int CONNECTIVITY::getNumberOfTypesWithPoly(MED_EN::medEntityMesh Entity) const
2477 {
2478   SCRUTE_MED(_entity);
2479   SCRUTE_MED(Entity);
2480   if (_entity==Entity)
2481     {
2482       SCRUTE_MED(_numberOfTypes);
2483       SCRUTE_MED(getNumberOfPolyType());
2484       return _numberOfTypes+getNumberOfPolyType();
2485     }
2486   else if (_constituent!=NULL)
2487     return _constituent->getNumberOfTypesWithPoly(Entity);
2488   else
2489     return 0;
2490 }
2491
2492 int CONNECTIVITY::getNumberOfPolyType()  const
2493 {
2494   if (_entity==MED_CELL && _entityDimension==3)
2495     {
2496       if(getNumberOfPolyhedron()>0)
2497         return 1;
2498     }
2499   else if ((_entity==MED_CELL && _entityDimension==2) || (_entity==MED_FACE))
2500     if(getNumberOfPolygons()>0)
2501       return 1;
2502   return 0;
2503 }
2504
2505 int CONNECTIVITY::getNumberOfElementOfPolyType(MED_EN::medEntityMesh Entity)  const
2506 {
2507   if(_entity==Entity)
2508     {
2509       if (_entity==MED_CELL && _entityDimension==3)
2510         return getNumberOfPolyhedron();
2511       else if ((_entity==MED_CELL && _entityDimension==2) || (_entity==MED_FACE))
2512         return getNumberOfPolygons();
2513       return 0;
2514     }
2515   else
2516     {
2517       if(_constituent!=NULL)
2518         return _constituent->getNumberOfElementOfPolyType(Entity);
2519       else
2520         //throw MEDEXCEPTION("_constituent required to evaluate getNumberOfElementOfPolyType");
2521         return 0;
2522     }
2523 }
2524
2525 /*
2526   Method equivalent to getGeometricTypes except that it includes not only classical Types but polygons/polyhedra also.
2527   WARNING the returned array MUST be deallocated.
2528  */
2529 MED_EN::medGeometryElement * CONNECTIVITY::getGeometricTypesWithPoly(MED_EN::medEntityMesh Entity) const throw (MEDEXCEPTION)
2530 {
2531   if(_entity==Entity)
2532     {
2533       int nbOfTypesTotal=getNumberOfTypesWithPoly(Entity);
2534       int nbOfTypesWithoutPoly=getNumberOfTypes(Entity);
2535       medGeometryElement *ret=new medGeometryElement[nbOfTypesTotal];
2536       memcpy(ret,getGeometricTypes(Entity),nbOfTypesWithoutPoly*sizeof(medGeometryElement));
2537       if(nbOfTypesTotal>nbOfTypesWithoutPoly)
2538         {
2539           if (Entity==MED_CELL && _entityDimension==3)
2540             ret[nbOfTypesWithoutPoly]=MED_EN::MED_POLYHEDRA;
2541           else if((Entity==MED_CELL && _entityDimension==2) || (Entity==MED_FACE))
2542             ret[nbOfTypesWithoutPoly]=MED_EN::MED_POLYGON;
2543         }
2544       return ret;
2545     }
2546   else
2547     {
2548       if(_constituent)
2549         return _constituent->getGeometricTypesWithPoly(Entity);
2550                         else
2551                                 return 0;
2552     }
2553 }
2554
2555 /*
2556   Method used in CalculateDescendingConnectivity. So it's typically a private method.
2557   The aim of this method is to hide to CalculateDescendingConnectivity algorithm the fact that in reverse connectivity polygons and polyhedrons
2558   are not in separate data structure, contrary to not reverse connectivities.
2559  */
2560 int CONNECTIVITY::getIndexOfEndClassicElementInReverseNodal(const int *reverseNodalValue, const int *reverseNodalIndex, int rk)  const
2561 {
2562   int nbOfLastElt=getNumberOf(_entity,MED_ALL_ELEMENTS);
2563   int ret=reverseNodalIndex[rk];
2564   for(int i=reverseNodalIndex[rk];i<reverseNodalIndex[rk+1];i++)
2565     {
2566       if(reverseNodalValue[i-1]<=nbOfLastElt)
2567         ret++;
2568     }
2569   return ret;
2570 }
2571
2572 /*
2573   Method that inverts the face with faceId 'faceId' in the data structure. If it's a polygon face 'faceId' is a value between 1 and nbOfPolygons.
2574   This method has to be applied on a instance of MEDMEM::CONNECTIVITY with _entityDimension==3.
2575   WARNING calculateDescendingConnectivity must have been called before.
2576  */
2577 void CONNECTIVITY::invertConnectivityForAFace(int faceId, const int *nodalConnForFace, bool polygonFace)
2578 {
2579   // we have always 2 neighbourings
2580   int cell1 = _reverseDescendingConnectivity->getIJ(faceId,1);
2581   int cell2 = _reverseDescendingConnectivity->getIJ(faceId,2);
2582
2583   if (cell2 != 0) { // we are not on border, update compulsary. If we aren't on border no update necessary so WARNING because user specified a bad oriented face
2584     // Updating _reverseDescendingConnectivity
2585     _reverseDescendingConnectivity->setIJ(faceId,1,cell2);
2586     _reverseDescendingConnectivity->setIJ(faceId,2,cell1);
2587     // Updating _constituent->_nodal because of reversity
2588     const int *descendingNodalIndex=(!polygonFace)?_constituent->_nodal->getIndex():_constituent->_polygonsNodal->getIndex();
2589     MEDSKYLINEARRAY *currentNodal=(!polygonFace)?_constituent->_nodal:_constituent->_polygonsNodal;
2590     int faceIdRelative=(!polygonFace)?faceId:faceId-getNumberOf(MED_FACE,MED_ALL_ELEMENTS);
2591     for(int iarray=1;iarray<=(descendingNodalIndex[faceIdRelative]-descendingNodalIndex[faceIdRelative-1]);iarray++)
2592       currentNodal->setIJ(faceIdRelative,iarray,nodalConnForFace[iarray-1]);
2593
2594     // Updating _descending for cell1 and cell2
2595     const int NB_OF_CELLS_SHARING_A_FACE=2;
2596     int cellsToUpdate[NB_OF_CELLS_SHARING_A_FACE]; cellsToUpdate[0]=cell1; cellsToUpdate[1]=cell2;
2597     for(int curCell=0;curCell<NB_OF_CELLS_SHARING_A_FACE;curCell++)
2598       {
2599         int cell=cellsToUpdate[curCell];
2600         bool polyhCell=(getElementTypeWithPoly(MED_CELL,cell)==MED_POLYHEDRA);
2601         if(polyhCell)
2602           cell-=getNumberOf(MED_CELL,MED_ALL_ELEMENTS);
2603         const int *newDescendingIndex=(!polyhCell)?_descending->getIndex():_polyhedronDescending->getIndex();
2604         MEDSKYLINEARRAY *currentDescending=(!polyhCell)?_descending:_polyhedronDescending;
2605         for(int iface=newDescendingIndex[cell-1];iface<newDescendingIndex[cell];iface++)
2606           {
2607             int curValue=currentDescending->getIndexValue(iface);
2608             if (abs(curValue)==faceId)
2609               currentDescending->setIndexValue(iface,-curValue);
2610           }
2611       }
2612   }
2613 }
2614
2615 /*
2616   Method with 2 output : the connectivity required and its length 'lgth'. This method gives the connectivity independently it is a polygons/polyhedrons or normal element.
2617  */
2618 const int * CONNECTIVITY::getConnectivityOfAnElementWithPoly(MED_EN::medConnectivity ConnectivityType, MED_EN::medEntityMesh Entity, int Number, int &lgth)
2619 {
2620   if(Entity==MED_EN::MED_NODE)
2621     throw  MEDEXCEPTION("No connectivity attached to a node entity");
2622   if(Entity==_entity)
2623     {
2624       if(_entity==MED_EDGE && _entityDimension==1)
2625         {
2626           const int * newConstituentValue = 0;
2627           const int * newConstituentIndex = 0;
2628           newConstituentValue = getConnectivity(ConnectivityType,Entity,MED_ALL_ELEMENTS);
2629           newConstituentIndex = getConnectivityIndex(ConnectivityType,Entity);
2630           lgth=newConstituentIndex[Number]-newConstituentIndex[Number-1];
2631           return newConstituentValue+newConstituentIndex[Number-1]-1;
2632         }
2633       int nbOfClassicalElements=getNumberOf(_entity,MED_ALL_ELEMENTS);
2634       if((_entity==MED_FACE) || (_entity==MED_CELL && _entityDimension==2))
2635         {
2636           const int * newConstituentValue = 0;
2637           const int * newConstituentIndex = 0;
2638           if(Number<=nbOfClassicalElements && nbOfClassicalElements!=0)
2639             {
2640               newConstituentValue = getConnectivity(ConnectivityType,Entity,MED_ALL_ELEMENTS);
2641               newConstituentIndex = getConnectivityIndex(ConnectivityType,Entity);
2642               lgth=newConstituentIndex[Number]-newConstituentIndex[Number-1];
2643               return newConstituentValue+newConstituentIndex[Number-1]-1;
2644             }
2645           else
2646             {
2647               int localNumber=Number-nbOfClassicalElements-1;
2648               if(localNumber<getNumberOfPolygons())
2649                 {
2650                   newConstituentValue = getPolygonsConnectivity(ConnectivityType,Entity);
2651                   newConstituentIndex = getPolygonsConnectivityIndex(ConnectivityType,Entity);
2652                   lgth=newConstituentIndex[localNumber+1]-newConstituentIndex[localNumber];
2653                   return newConstituentValue+newConstituentIndex[localNumber]-1;
2654                 }
2655               else
2656                 throw  MEDEXCEPTION("Unknown number");
2657             }
2658         }
2659       else if(_entity==MED_CELL && _entityDimension==3)
2660         {
2661           const int * newConstituentValue = 0;
2662           const int * newConstituentIndex = 0;
2663           if(Number<=nbOfClassicalElements)
2664             {
2665               newConstituentValue = getConnectivity(ConnectivityType,Entity,MED_ALL_ELEMENTS);
2666               newConstituentIndex = getConnectivityIndex(ConnectivityType,Entity);
2667               lgth=newConstituentIndex[Number]-newConstituentIndex[Number-1];
2668               return newConstituentValue+newConstituentIndex[Number-1]-1;
2669             }
2670           else
2671             {
2672               int localNumber=Number-nbOfClassicalElements-1;
2673               if(localNumber<getNumberOfPolyhedron())
2674                 {
2675                   if(ConnectivityType==MED_NODAL)
2676                     throw  MEDEXCEPTION("NODAL Connectivity for a polyhedron not directly accessible.\n Use getPolyhedronNodal and getPolyhedronFaces instead");
2677 //                newConstituentValue = _polyhedronDescending->getValue();
2678 //                newConstituentIndex = _polyhedronDescending->getIndex();
2679                   newConstituentValue = getPolyhedronConnectivity( ConnectivityType );
2680                   newConstituentIndex = getPolyhedronIndex( ConnectivityType );
2681                   lgth=newConstituentIndex[localNumber+1]-newConstituentIndex[localNumber];
2682                   return newConstituentValue+newConstituentIndex[localNumber]-1;
2683                 }
2684               else
2685                 throw  MEDEXCEPTION("Unknown number");
2686             }
2687         }
2688       else
2689         throw MEDEXCEPTION("No connectivity available");
2690     }
2691   else
2692     {
2693       if(_constituent==NULL)
2694         calculateDescendingConnectivity();
2695       return _constituent->getConnectivityOfAnElementWithPoly(ConnectivityType,Entity,Number,lgth);
2696     }
2697 }
2698
2699 int CONNECTIVITY::getNumberOfElementsWithPoly(MED_EN::medEntityMesh Entity, MED_EN::medGeometryElement Type) const
2700 {
2701   if(Entity==MED_EN::MED_NODE)
2702     return _numberOfNodes;
2703   if(Entity==_entity)
2704     {
2705       if (Type == MED_ALL_ELEMENTS)
2706         return  getNumberOfElementOfPolyType(_entity)+ getNumberOf(_entity,Type) ;
2707
2708       if(Type==MED_POLYGON || Type==MED_POLYHEDRA)
2709         {
2710           if(Type==MED_POLYGON && Entity==MED_CELL && _entityDimension==3 || Type==MED_POLYHEDRA && Entity==MED_FACE)
2711             return 0;
2712           return getNumberOfElementOfPolyType(_entity);
2713         }
2714       else
2715         return getNumberOf(_entity,Type);
2716     }
2717   else
2718     {
2719       if(_constituent)
2720         return _constituent->getNumberOfElementsWithPoly(Entity,Type);
2721       else
2722         //throw MEDEXCEPTION("CONNECTIVITY::getNumberOfElementsWithPoly : _constituent needed");
2723         return 0;
2724     }
2725 }
2726
2727 /*!
2728   Perform a deep comparison of the 2 connectivities in NODAL mode and on all elements.
2729 */
2730 bool CONNECTIVITY::deepCompare(const CONNECTIVITY& other) const
2731 {
2732   CONNECTIVITY* temp=(CONNECTIVITY* )this;
2733   const int *conn1=temp->getConnectivity(MED_NODAL,_entity,MED_ALL_ELEMENTS);
2734   int size1=temp->getConnectivityLength(MED_NODAL,_entity,MED_ALL_ELEMENTS);
2735   temp=(CONNECTIVITY* )(&other);
2736   const int *conn2=temp->getConnectivity(MED_NODAL,_entity,MED_ALL_ELEMENTS);
2737   int size2=temp->getConnectivityLength(MED_NODAL,_entity,MED_ALL_ELEMENTS);
2738   if(size1!=size2)
2739     return false;
2740   bool ret=true;
2741   for(int i=0;i<size1 && ret;i++)
2742     {
2743       ret=(conn1[i]==conn2[i]);
2744     }
2745   return ret;
2746 }