src/MEDMEM/MEDMEM_Support.cxx

   1 /*
   2  File Support.cxx
   3  $Header$
   4 */
   5
   6 #include <set>
   7 #include <algorithm>
   8 #include <list>
   9
  10 #include "MEDMEM_DriversDef.hxx"
  11 #include "MEDMEM_Support.hxx"
  12 #include "MEDMEM_Mesh.hxx"
  13
  14 using namespace std;
  15 using namespace MED_EN;
  16 using namespace MEDMEM;
  17
  18 #define MED_NBR_GEOMETRIE_MAILLE 15
  19
  20 /* This class is a generic class for family and group */
  21
  22 /*!
  23   Constructor.
  24 */
  25 //--------------------------------------------------------------------------
  26 SUPPORT::SUPPORT(): _name(""),  _description("None"), _mesh((MESH*)NULL),
  27                     _entity(MED_CELL), _numberOfGeometricType(0),
  28                     _geometricType((medGeometryElement*)NULL),
  29                     _numberOfGaussPoint((int*)NULL),
  30                     _isOnAllElts(false),
  31                     _numberOfElements((int*)NULL),
  32                     _totalNumberOfElements(0),
  33                     _number((MEDSKYLINEARRAY*)NULL)
  34 //--------------------------------------------------------------------------
  35 {
  36     MESSAGE("SUPPORT::SUPPORT()");
  37 };
  38
  39 /*!
  40   Constructor.
  41 */
  42 //--------------------------------------------------------------------------
  43 SUPPORT::SUPPORT(MESH* Mesh, string Name/*=""*/, medEntityMesh Entity/*=MED_CELL*/):
  44                 _name(Name), _description("None"), _mesh(Mesh), _entity(Entity),
  45                 _numberOfGeometricType(0),
  46                 _geometricType((medGeometryElement*)NULL),
  47                 _numberOfGaussPoint((int*)NULL),
  48                 _isOnAllElts(true),
  49                 _numberOfElements((int*)NULL),
  50                 _totalNumberOfElements(0),
  51                 _number((MEDSKYLINEARRAY*)NULL)
  52 //--------------------------------------------------------------------------
  53 {
  54   MESSAGE("SUPPORT::SUPPORT(MESH*Mesh,string Name,medEntityMesh Entity)");
  55   update() ;
  56 };
  57
  58 /*!
  59   Copy constructor.
  60 */
  61 //--------------------------------------------------------------------------
  62 SUPPORT::SUPPORT(const SUPPORT & m)
  63 //--------------------------------------------------------------------------
  64 {
  65   const char * LOC = "SUPPORT::SUPPORT(SUPPORT & m) : " ;
  66   BEGIN_OF(LOC) ;
  67
  68   _name = m._name;
  69   _description = m._description;
  70   _mesh = m._mesh; // on recopie uniquement l'adresse
  71   _entity = m._entity;
  72   _numberOfGeometricType = m._numberOfGeometricType;
  73   if (m._geometricType != NULL)
  74     {
  75       _geometricType = new medGeometryElement[m._numberOfGeometricType];
  76       memcpy(_geometricType,m._geometricType,m._numberOfGeometricType*sizeof(medGeometryElement));
  77     }
  78   else
  79     _geometricType = (medGeometryElement *) NULL;
  80   if (m._numberOfGaussPoint != NULL)
  81     {
  82       _numberOfGaussPoint = new int[m._numberOfGeometricType];
  83       memcpy(_numberOfGaussPoint,m._numberOfGaussPoint,m._numberOfGeometricType*sizeof(int));
  84     }
  85   else
  86     _numberOfGaussPoint = (int *) NULL;
  87   _isOnAllElts = m._isOnAllElts;
  88   if (m._numberOfElements != NULL)
  89     {
  90       _numberOfElements = new int[_numberOfGeometricType];
  91       memcpy(_numberOfElements,m._numberOfElements,_numberOfGeometricType*sizeof(int));
  92     }
  93   else
  94     _numberOfElements = (int *) NULL;
  95   _totalNumberOfElements = m._totalNumberOfElements;
  96   if (m._isOnAllElts == false)
  97     _number = new MEDSKYLINEARRAY(* m._number);
  98   else
  99     _number = (MEDSKYLINEARRAY *) NULL;
 100
 101   END_OF(LOC) ;
 102 };
 103
 104
 105 /*!
 106   Destructor.
 107 */
 108 //-----------------
 109 SUPPORT::~SUPPORT()
 110 //-----------------
 111 {
 112   MESSAGE("Destructeur ~SUPPORT()");
 113   clearDataOnNumbers();
 114 }
 115
 116 /*!
 117   operator <<.
 118 */
 119 //--------------------------------------------------
 120 ostream & MEDMEM::operator<<(ostream &os, const SUPPORT &my)
 121 //--------------------------------------------------
 122 {
 123   os << "Name : "<< my.getName() << endl ;
 124   os << "Description : "<< my.getDescription() << endl ;
 125   os << "Mesh name : ";
 126   if (my.getMesh() == NULL)
 127     os << "ERROR : Mesh not defined !" << endl ;
 128   else
 129     os << my._mesh->getName() << endl ;
 130   os << "Entity : "<< my._entity << endl;
 131   os << "Entity list : "<< endl;
 132   if (!(my._isOnAllElts)) {
 133     int numberoftypes = my._numberOfGeometricType ;
 134     os << "NumberOfTypes : "<<numberoftypes<<endl;
 135     medGeometryElement * types = my._geometricType;
 136     for (int j=0;j<numberoftypes;j++) {
 137       int numberOfElements = my._numberOfElements[j];
 138       os << "    * Type "<<types[j]<<" : there is(are) "<<numberOfElements<<" element(s) :" << endl;
 139 //       const int * number = my.getNumber(types[j]);
 140 //       SCRUTE(number);
 141 //       os << " --> ";
 142 //       for (int k=0; k<numberOfElements;k++)
 143 //      os << number[k] << " ";
 144 //       os << endl ;
 145     }
 146   } else
 147     os << "Is on all entities !"<< endl;
 148
 149   return os ;
 150 }
 151
 152 /*!
 153   Updade the SUPPORT attributs with rigth MESH information.
 154
 155   It has an effect only if SUPPORT is on all elements.
 156
 157   No more need in future release.
 158 */
 159 //-------------------
 160 void SUPPORT::update()
 161 //-------------------
 162 {
 163   const char * LOC = "SUPPORT::update() : " ;
 164   BEGIN_OF(LOC) ;
 165
 166   if (_isOnAllElts) {
 167     if (_entity == MED_NODE) {
 168       _numberOfGeometricType=1 ;
 169       _geometricType=new medGeometryElement[1] ;
 170       _geometricType[0]=MED_NONE ;
 171       _numberOfElements = new int[1] ;
 172       _numberOfElements[0]=_mesh->getNumberOfNodes();
 173       _totalNumberOfElements=_numberOfElements[0];
 174       _numberOfGaussPoint = new int[1] ;
 175       _numberOfGaussPoint[0]=1;
 176     } else { // we duplicate information from _mesh
 177       _numberOfGeometricType=_mesh->getNumberOfTypesWithPoly(_entity);
 178       if (_geometricType == (medGeometryElement *) NULL)
 179         _geometricType=_mesh->getTypesWithPoly(_entity);
 180       else
 181         memcpy(_geometricType,_mesh->getTypes(_entity),_numberOfGeometricType*sizeof(medGeometryElement));
 182       if (_numberOfElements == (int *) NULL)
 183         _numberOfElements = new int[_numberOfGeometricType] ;
 184       if (_numberOfGaussPoint == (int *) NULL)
 185         _numberOfGaussPoint = new int[_numberOfGeometricType] ;
 186       _totalNumberOfElements=0;
 187       for (int i=0;i<_numberOfGeometricType;i++) {
 188         _numberOfElements[i]=_mesh->getNumberOfElementsWithPoly(_entity,_geometricType[i]) ;
 189         _totalNumberOfElements+=_numberOfElements[i];
 190         _numberOfGaussPoint[i]=1 ;
 191       }
 192     }
 193   }
 194   END_OF(LOC);
 195 };
 196
 197 /*!
 198   Blend the given SUPPORT mySupport into the calling object SUPPORT.
 199 */
 200 //-------------------
 201 void SUPPORT::blending(SUPPORT * mySupport) throw (MEDEXCEPTION)
 202 //-------------------
 203 {
 204   const char * LOC="SUPPORT::blending(SUPPORT *) : ";
 205   BEGIN_OF(LOC);
 206   if (_entity!=mySupport->getEntity())
 207     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Entities are different !"));
 208   if(!(*_mesh == *mySupport->getMesh()))
 209     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Mesh are different !"));
 210   if(_isOnAllElts)
 211     return;
 212   if(mySupport->isOnAllElements())
 213     {
 214       *this=*mySupport;
 215       return;
 216     }
 217   if(mySupport->_totalNumberOfElements==0)
 218     return;
 219   const int *ids=getNumber(MED_ALL_ELEMENTS);
 220   set<int> idsSet(ids,ids+getNumberOfElements(MED_ALL_ELEMENTS));
 221   const int *idsMySupport=mySupport->getNumber(MED_ALL_ELEMENTS);
 222   int mySupportSize=mySupport->getNumberOfElements(MED_ALL_ELEMENTS);
 223   set<int>::iterator iter;
 224   for(int i=0;i<mySupportSize;i++)
 225     idsSet.insert(idsMySupport[i]);
 226   int size=idsSet.size();
 227   if(size!=0)
 228     {
 229       list<int> idsList;
 230       for(iter=idsSet.begin();iter!=idsSet.end();iter++)
 231         idsList.push_back(*iter);
 232       if(_entity==MED_NODE)
 233         fillFromNodeList(idsList);
 234       else
 235         fillFromElementList(idsList);
 236     }
 237   else
 238     clearDataOnNumbers();
 239   END_OF(LOC);
 240 }
 241
 242 /*!
 243     This function allows the user to set a support not on all entities Entity,
 244     it should be used after an initialisation with the constructor
 245     SUPPORT(MESH* Mesh, string Name="", medEntityMesh Entity=MED_CELL) and
 246     after the call to the function setAll(false).
 247     It allocates and initialises all the attributs of the class SUPPORT.
 248  */
 249
 250 //-------------------
 251 void SUPPORT::setpartial(string Description, int NumberOfGeometricType,
 252                          int TotalNumberOfElements,
 253                          medGeometryElement *GeometricType,
 254                          int *NumberOfElements, int *NumberValue)
 255 //-------------------
 256 {
 257   const char * LOC = "SUPPORT::setpartial(string , int , int , medGeometryElement * , int * , int *) : " ;
 258   BEGIN_OF(LOC) ;
 259
 260   _isOnAllElts = false ;
 261
 262   _description=Description;
 263
 264   _numberOfGeometricType=NumberOfGeometricType;
 265
 266   if (_geometricType!=NULL) delete[] _geometricType ;
 267   _geometricType = new medGeometryElement[NumberOfGeometricType];
 268   if (_numberOfElements!=NULL) delete[] _numberOfElements ;
 269   _numberOfElements = new int[NumberOfGeometricType];
 270   _totalNumberOfElements = TotalNumberOfElements;
 271   if (_numberOfGaussPoint!=NULL) delete[] _numberOfGaussPoint ;
 272   _numberOfGaussPoint = new int[NumberOfGeometricType];
 273   int * index = new int[_numberOfGeometricType+1];
 274   index[0]=1;
 275   int elemDim = -1;
 276   for (int i=0;i<_numberOfGeometricType;i++) {
 277     if(GeometricType[i]/100 != elemDim)
 278       if(i==0)
 279         elemDim=GeometricType[i]/100;
 280       else
 281         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"unhomogeneous geometric types (dimension) !"));
 282     _geometricType[i] = GeometricType[i] ;
 283     _numberOfElements[i] = NumberOfElements[i] ;
 284     _numberOfGaussPoint[i] = 1 ;
 285     index[i+1] = index[i]+NumberOfElements[i] ;
 286   }
 287
 288   if (_number!=NULL) delete _number ;
 289   _number = new MEDSKYLINEARRAY(_numberOfGeometricType,_totalNumberOfElements,index,NumberValue);
 290
 291   delete[] index ;
 292
 293   END_OF(LOC);
 294 };
 295
 296
 297 /*!
 298   This method gets the boundary elements of the mesh. The support has to be
 299   build using the constructor SUPPORT(MESH *,string, medEntityMesh) or
 300   SUPPORT() followed by setMesh(MESH*) setName(string) and
 301   setEntity(medEntityMesh) before using this method.
 302 */
 303 //-------------------
 304 void SUPPORT::getBoundaryElements() throw (MEDEXCEPTION)
 305 //-------------------
 306 {
 307   const char * LOC = "SUPPORT::getBoundaryElements() : " ;
 308   BEGIN_OF(LOC) ;
 309
 310   if (_mesh == (MESH*)NULL) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"You shlould use the appropriate SUPPORT Constructor before calling this method"));
 311
 312   int spaceDimension = _mesh->getSpaceDimension();
 313
 314   if (spaceDimension == 3)
 315     if (_entity != MED_FACE)
 316       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<_entity<<" !"));
 317   if (spaceDimension == 2)
 318     if (_entity != MED_EDGE)
 319       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<_entity<<" !"));
 320
 321   setAll(false);
 322
 323   const int * myConnectivityValue = _mesh->getReverseConnectivity(MED_DESCENDING) ;
 324   const int * myConnectivityIndex = _mesh->getReverseConnectivityIndex(MED_DESCENDING) ;
 325   int numberOf = _mesh->getNumberOfElements(_entity,MED_ALL_ELEMENTS) ;
 326   list<int> myElementsList ;
 327   int size = 0 ;
 328   SCRUTE(numberOf) ;
 329   for (int i=0 ; i<numberOf; i++)
 330     if (myConnectivityValue[myConnectivityIndex[i]] == 0) {
 331       SCRUTE(i+1) ;
 332       myElementsList.push_back(i+1) ;
 333       size++ ;
 334     }
 335   SCRUTE(size) ;
 336   // Well, we must know how many geometric type we have found
 337   int * myListArray = new int[size] ;
 338   int id = 0 ;
 339   list<int>::iterator myElementsListIt ;
 340   for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
 341     myListArray[id]=(*myElementsListIt) ;
 342     SCRUTE(id);
 343     SCRUTE(myListArray[id]);
 344     id ++ ;
 345   }
 346
 347   int numberOfGeometricType ;
 348   medGeometryElement* geometricType ;
 349   int * numberOfGaussPoint ;
 350   int * geometricTypeNumber ;
 351   int * numberOfElements ;
 352   //MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY() ;
 353   int * mySkyLineArrayIndex ;
 354
 355   int numberOfType = _mesh->getNumberOfTypes(_entity) ;
 356   if (numberOfType == 1) { // wonderfull : it's easy !
 357     numberOfGeometricType = 1 ;
 358     geometricType = new medGeometryElement[1] ;
 359     const medGeometryElement *  allType = _mesh->getTypes(_entity);
 360     geometricType[0] = allType[0] ;
 361     numberOfGaussPoint = new int[1] ;
 362     numberOfGaussPoint[0] = 1 ;
 363     geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
 364     geometricTypeNumber[0] = 0 ;
 365     numberOfElements = new int[1] ;
 366     numberOfElements[0] = size ;
 367     mySkyLineArrayIndex = new int[2] ;
 368     mySkyLineArrayIndex[0]=1 ;
 369     mySkyLineArrayIndex[1]=1+size ;
 370   }
 371   else {// hemmm
 372     map<medGeometryElement,int> theType ;
 373     for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
 374       medGeometryElement myType = _mesh->getElementType(_entity,*myElementsListIt) ;
 375       if (theType.find(myType) != theType.end() )
 376         theType[myType]+=1 ;
 377       else
 378         theType[myType]=1 ;
 379     }
 380     numberOfGeometricType = theType.size() ;
 381     geometricType = new medGeometryElement[numberOfGeometricType] ;
 382     //const medGeometryElement *  allType = _mesh->getTypes(_entity); !! UNUSED VARIABLE !!
 383     numberOfGaussPoint = new int[numberOfGeometricType] ;
 384     geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
 385     numberOfElements = new int[numberOfGeometricType] ;
 386     mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
 387     int index = 0 ;
 388     mySkyLineArrayIndex[0]=1 ;
 389     map<medGeometryElement,int>::iterator theTypeIt ;
 390     for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
 391       geometricType[index] = (*theTypeIt).first ;
 392       numberOfGaussPoint[index] = 1 ;
 393       geometricTypeNumber[index] = 0 ;
 394       numberOfElements[index] = (*theTypeIt).second ;
 395       mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfElements[index] ;
 396       index++ ;
 397     }
 398   }
 399   //mySkyLineArray->setMEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
 400   MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
 401
 402   setNumberOfGeometricType(numberOfGeometricType) ;
 403   //  setGeometricType(geometricType) ;
 404   //  setNumberOfGaussPoint(numberOfGaussPoint) ;
 405   for (int i=0;i<numberOfGeometricType;i++)
 406     {
 407       _numberOfGaussPoint[i] = numberOfGaussPoint[i];
 408       _geometricType[i] = geometricType[i];
 409     }
 410
 411   setNumberOfElements(numberOfElements) ;
 412   setTotalNumberOfElements(size) ;
 413   //  setNumber(mySkyLineArray) ;
 414
 415   _number = new MEDSKYLINEARRAY(numberOfGeometricType,size);
 416
 417   _number->setIndex(mySkyLineArrayIndex);
 418
 419   for (int i=0;i<size;i++)
 420     {
 421       _number->setIndexValue(i+1,myListArray[i]);
 422     }
 423
 424   delete[] numberOfElements;
 425   delete[] geometricTypeNumber;
 426   delete[] numberOfGaussPoint;
 427   delete[] geometricType;
 428   delete[] mySkyLineArrayIndex;
 429   delete[] myListArray;
 430   delete mySkyLineArray;
 431
 432   END_OF(LOC) ;
 433 }
 434
 435 /*!
 436   intersect the given SUPPORT mySupport into the calling SUPPORT object.
 437 */
 438 //-------------------
 439 void SUPPORT::intersecting(SUPPORT * mySupport) throw (MEDEXCEPTION)
 440 {
 441   const char * LOC="SUPPORT::intersecting(SUPPORT *) : ";
 442   BEGIN_OF(LOC);
 443   if (_entity!=mySupport->getEntity())
 444     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Entities are different !"));
 445   if(!(*_mesh == *mySupport->getMesh()))
 446     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Mesh are different !"));
 447   if(mySupport->isOnAllElements())
 448     return;
 449   if(_isOnAllElts)
 450     {
 451       *this=*mySupport;
 452       return;
 453     }
 454   if(_totalNumberOfElements==0)
 455     return;
 456   const int *ids=getNumber(MED_ALL_ELEMENTS);
 457   set<int> idsSet(ids,ids+getNumberOfElements(MED_ALL_ELEMENTS));
 458   const int *idsMySupport=mySupport->getNumber(MED_ALL_ELEMENTS);
 459   int mySupportSize=mySupport->getNumberOfElements(MED_ALL_ELEMENTS);
 460   set<int> idsSetMySupport(idsMySupport,idsMySupport+mySupportSize);
 461   set<int>::iterator iter;
 462   list<int> idsList;
 463   for(iter=idsSet.begin();iter!=idsSet.end();iter++)
 464     if(idsSetMySupport.find(*iter)!=idsSetMySupport.end())
 465       idsList.push_back(*iter);
 466   int size=idsSet.size();
 467   if(size!=0)
 468     {
 469       if(_entity==MED_NODE)
 470         fillFromNodeList(idsList);
 471       else
 472         fillFromElementList(idsList);
 473     }
 474   else
 475     {
 476       clearDataOnNumbers();
 477     }
 478   END_OF(LOC);
 479 };
 480
 481 /*!
 482   operator = perform et deep copy except for attribute _mesh
 483  */
 484 SUPPORT& MEDMEM::SUPPORT::operator=(const SUPPORT &other)
 485 {
 486   const char * LOC = "SUPPORT::operator=(const SUPPORT & m) : ";
 487   BEGIN_OF(LOC);
 488
 489   _name=other._name;
 490   _description=other._description;
 491   _mesh=other._mesh;
 492   _entity=other._entity;
 493   _isOnAllElts=other._isOnAllElts;
 494   clearDataOnNumbers();
 495   _numberOfGeometricType=other._numberOfGeometricType;
 496   _totalNumberOfElements=other._totalNumberOfElements;
 497
 498   if(other._geometricType != NULL)
 499     {
 500       _geometricType=new medGeometryElement[other._numberOfGeometricType];
 501       memcpy(_geometricType,other._geometricType,other._numberOfGeometricType*sizeof(medGeometryElement));
 502     }
 503   if (other._numberOfGaussPoint != NULL)
 504     {
 505       _numberOfGaussPoint=new int[other._numberOfGeometricType];
 506       memcpy(_numberOfGaussPoint,other._numberOfGaussPoint,other._numberOfGeometricType*sizeof(int));
 507     }
 508   if(other._numberOfElements != NULL)
 509     {
 510       _numberOfElements=new int[_numberOfGeometricType];
 511       memcpy(_numberOfElements,other._numberOfElements,_numberOfGeometricType*sizeof(int));
 512     }
 513   if(!_isOnAllElts)
 514     {
 515       _number=new MEDSKYLINEARRAY(* other._number);
 516     }
 517   return *this;
 518 }
 519
 520 /*!
 521   Method that cleans up all the fields related to _numbers. Defined for code factorization.
 522  */
 523 //--------------------------------------------------
 524 void MEDMEM::SUPPORT::clearDataOnNumbers()
 525 //--------------------------------------------------
 526 {
 527   _numberOfGeometricType=0;
 528   _totalNumberOfElements=0;
 529   if(_geometricType)
 530     {
 531       delete [] _geometricType;
 532       _geometricType=(medGeometryElement *) NULL;
 533     }
 534   if(_numberOfGaussPoint)
 535     {
 536       delete [] _numberOfGaussPoint;
 537       _numberOfGaussPoint=(int *) NULL;
 538     }
 539   if(_numberOfElements)
 540     {
 541       delete [] _numberOfElements;
 542       _numberOfElements=(int *) NULL;
 543     }
 544   if(_number)
 545     {
 546       delete _number;
 547       _number=(MEDSKYLINEARRAY *) NULL;
 548     }
 549 }
 550
 551 /*!
 552   operator == This operator does not compare attributs _name and _description.
 553 */
 554 //--------------------------------------------------
 555 bool MEDMEM::SUPPORT::operator == (const SUPPORT &support) const
 556 //--------------------------------------------------
 557 {
 558   const char * LOC = "bool SUPPORT::operator ==(const SUPPORT &support) const : ";
 559
 560   BEGIN_OF(LOC);
 561
 562   bool operatorReturn = false;
 563
 564   operatorReturn = (*_mesh == *support._mesh) && (_entity == support._entity) &&
 565     (_numberOfGeometricType == support._numberOfGeometricType) &&
 566     ((_isOnAllElts && support._isOnAllElts) || (!_isOnAllElts && !support._isOnAllElts)) &&
 567     (_totalNumberOfElements == support._totalNumberOfElements);
 568
 569   if (operatorReturn)
 570     {
 571       if (!_isOnAllElts)
 572         {
 573           for (int i=0; i<_numberOfGeometricType; i++)
 574             {
 575               operatorReturn = operatorReturn &&
 576                 (_geometricType[i] == support._geometricType[i]) &&
 577                 (_numberOfElements[i] == support._numberOfElements[i]) &&
 578                 (_numberOfGaussPoint[i] == support._numberOfGaussPoint[i]);
 579
 580               if (operatorReturn)
 581                 {
 582                   for (int j=0; j<_numberOfElements[i]; j++)
 583                     {
 584                       operatorReturn = operatorReturn &&
 585                         (getNumber(_geometricType[i])[j] ==
 586                          support.getNumber(_geometricType[i])[j]);
 587                     }
 588                 }
 589             }
 590         }
 591     }
 592
 593   END_OF(LOC);
 594
 595   return operatorReturn;
 596 };
 597
 598 void SUPPORT::changeElementsNbs(MED_EN::medEntityMesh entity, const int *renumberingFromOldToNew, int limitNbClassicPoly, const int *renumberingFromOldToNewPoly)
 599 {
 600   if(entity != _entity)
 601     throw MEDEXCEPTION("SUPPORT::changeElementsNbs : Renumbering on a mismatch entity");
 602   list<int> newNbs;
 603   if(!_isOnAllElts)
 604     {
 605       const int *oldNbs=_number->getValue();
 606       for(int i=0;i<_totalNumberOfElements;i++)
 607         {
 608           int globNb=oldNbs[i];
 609           if(globNb<=limitNbClassicPoly)
 610             newNbs.push_back(renumberingFromOldToNew[globNb-1]);
 611           else
 612             newNbs.push_back(renumberingFromOldToNewPoly[globNb-limitNbClassicPoly-1]);
 613         }
 614       newNbs.sort();
 615       fillFromElementList(newNbs);
 616     }
 617   else
 618     update();
 619 }
 620
 621 /*!
 622   operator == + in case false a test if coordinates and connectivity of _mesh and support->_mesh are the same
 623 */
 624 bool MEDMEM::SUPPORT::deepCompare(const SUPPORT &support) const
 625 {
 626   bool operatorReturn =(_entity == support._entity) &&
 627     (_numberOfGeometricType == support._numberOfGeometricType) &&
 628     ( (_isOnAllElts && support._isOnAllElts) || (!_isOnAllElts  && !support._isOnAllElts) ) &&
 629     (_totalNumberOfElements == support._totalNumberOfElements);
 630   if (operatorReturn)
 631     {
 632       if (!_isOnAllElts)
 633         {
 634           for (int i=0; i<_numberOfGeometricType && operatorReturn; i++)
 635             {
 636               operatorReturn = (_geometricType[i] == support._geometricType[i]) &&
 637                 (_numberOfElements[i] == support._numberOfElements[i]) &&
 638                 (_numberOfGaussPoint[i] == support._numberOfGaussPoint[i]);
 639               if (operatorReturn)
 640                 {
 641                   for (int j=0; j<_numberOfElements[i]; j++)
 642                     {
 643                       operatorReturn = (getNumber(_geometricType[i])[j] ==
 644                          support.getNumber(_geometricType[i])[j]);
 645                     }
 646                 }
 647             }
 648         }
 649     }
 650   if(operatorReturn)
 651     {
 652       if(!(*_mesh == *support._mesh))
 653         {
 654           return _mesh->deepCompare(*support._mesh);
 655         }
 656     }
 657   return operatorReturn;
 658 }
 659
 660 /*!
 661  States if this is included in other.
 662  */
 663 bool MEDMEM::SUPPORT::belongsTo(const SUPPORT& other, bool deepCompare) const
 664 {
 665   if(!(*_mesh == *other._mesh))
 666     {
 667       if(!deepCompare)
 668         return false;
 669       if(!_mesh->deepCompare(*other._mesh))
 670         return false;
 671     }
 672   if(_entity!=other._entity)
 673     return false;
 674   if(other._isOnAllElts)
 675     return true;
 676   if(_isOnAllElts && !other._isOnAllElts)
 677     return false;
 678   if(_numberOfGeometricType>other._numberOfGeometricType)
 679     return false;
 680   for(int i=0; i<_numberOfGeometricType; i++)
 681     {
 682       MED_EN::medGeometryElement curGeomType=_geometricType[i];
 683       int iOther=-1;
 684       for(int j=0; j<other._numberOfGeometricType; j++)
 685           if(other._geometricType[j]==curGeomType)
 686               iOther=j;
 687       if(iOther==-1)
 688         return false;
 689       if(_numberOfElements[i]>other._numberOfElements[iOther])
 690         return false;
 691       const int *numbers1=_number->getI(i+1);
 692       const int *numbers2=other._number->getI(iOther+1);
 693       for (int k=0; k<_numberOfElements[i]; k++)
 694         {
 695           bool found=false;
 696           for(int l=0;l<other._numberOfElements[iOther] && !found;l++)
 697             {
 698               if(numbers1[k]==numbers2[l])
 699                 found=true;
 700             }
 701           if(!found)
 702             return false;
 703         }
 704     }
 705   return true;
 706 }
 707 /*!
 708   Method used to sort array of id.
 709  */
 710 int compareId(const void *x, const void *y)
 711 {
 712   const int *x1=(const int *)x;
 713   const int *y1=(const int *)y;
 714   if(*x1<*y1)
 715     return -1;
 716   else if(*x1>*y1)
 717     return 1;
 718   else
 719     return 0;
 720 }
 721
 722 /*!
 723   performs a common operation : Sub builds a sorted int array that is obtained by supression of all ids contained
 724   in array defined by (idsToSuppress,lgthIdsToSuppress) from array [start ... end]
 725   Example sub(0,7,{1,2,5},3) => {0,3,4,6,7} - WARNING returned list should be deallocated !
 726  */
 727 list<int> *MEDMEM::SUPPORT::sub(int start,int end,const int *idsToSuppress,int lgthIdsToSuppress)
 728 {
 729   int size=end-start+1;
 730   int sizeRet=size-lgthIdsToSuppress;
 731   list<int> *ret;
 732   if(sizeRet<0)
 733     throw MEDEXCEPTION("MEDMEM::SUPPORT::sub");
 734   else if(sizeRet==0)
 735     {
 736       return 0;
 737     }
 738   if(idsToSuppress==0)
 739     {
 740       ret=new list<int>;
 741       for(int l=0;l<size;l++)
 742         ret->push_back(start+l);
 743       return ret;
 744     }
 745   ret=new list<int>;
 746   int *temp=new int[lgthIdsToSuppress];
 747   memcpy(temp,idsToSuppress,sizeof(int)*lgthIdsToSuppress);
 748   qsort(temp,lgthIdsToSuppress,sizeof(int),compareId);
 749   int k=0;
 750   for(int i=start;i<=end;i++)
 751     if(temp[k]!=i)
 752       ret->push_back(i);
 753     else
 754       k++;
 755   delete [] temp;
 756   return ret;
 757 }
 758
 759 /*!
 760   performs a common operation : Sub builds a sorted int array that is obtained by supression of all ids contained
 761   in array defined by (idsToSuppress,lgthIdsToSuppress) from array [start ... end]
 762   Example sub({1,3,4,5,6,7,9},7,{1,2,5},3) => {3,4,6,7,9}  - WARNING returned list should be deallocated !
 763  */
 764 list<int> *MEDMEM::SUPPORT::sub(const int *ids,int lgthIds,const int *idsToSuppress,int lgthIdsToSuppress)
 765 {
 766   list<int> *ret;
 767   int i,j=0;
 768   if(lgthIds<0)
 769     throw MEDEXCEPTION("MEDMEM::SUPPORT::sub");
 770   else if(lgthIds==0)
 771     return 0;
 772   ret=new list<int>;
 773   int *temp1=new int[lgthIds];
 774   memcpy(temp1,ids,sizeof(int)*lgthIds);
 775   qsort(temp1,lgthIds,sizeof(int),compareId);
 776   int *temp2=new int[lgthIdsToSuppress];
 777   memcpy(temp2,idsToSuppress,sizeof(int)*lgthIdsToSuppress);
 778   qsort(temp2,lgthIdsToSuppress,sizeof(int),compareId);
 779   for(i=0;i<lgthIds;)
 780     {
 781       if(j>=lgthIdsToSuppress)
 782           ret->push_back(temp1[i++]);
 783       else if(temp1[i]>temp2[j])
 784         j++;
 785       else if(temp1[i]<temp2[j])
 786         ret->push_back(temp1[i++]);
 787       else
 788         i++;
 789     }
 790   delete [] temp1;
 791   delete [] temp2;
 792   return ret;
 793 }
 794
 795 /*!
 796   returns a new SUPPORT (responsability to caller to destroy it)
 797   that is the complement to "this", lying on the same entity than "this".
 798  */
 799 SUPPORT *MEDMEM::SUPPORT::getComplement() const
 800 {
 801   SUPPORT *ret;
 802   const int nbOfElt=_mesh->getNumberOfElements(_entity,MED_EN::MED_ALL_ELEMENTS);
 803   int nbOfEltInSupp=getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
 804   if(_isOnAllElts || nbOfElt==nbOfEltInSupp)
 805     {
 806       ret=new SUPPORT;
 807       ret->setMesh(_mesh);
 808       ret->setEntity(_entity);
 809       string name="Complement of ";
 810       name+=_name;
 811       ret->setName(name);
 812       return ret;
 813     }
 814   const int *nbs=_number->getValue();
 815   list<int> *ids=sub(1,nbOfElt,nbs,nbOfEltInSupp);
 816   if(_entity==MED_EN::MED_NODE)
 817     ret=_mesh->buildSupportOnNodeFromElementList(*ids,_entity);
 818   else
 819     ret=_mesh->buildSupportOnElementsFromElementList(*ids,_entity);
 820   delete ids;
 821   return ret;
 822 }
 823
 824 /*!
 825   returns a new support the user should delete.
 826  */
 827 SUPPORT *MEDMEM::SUPPORT::substract(const SUPPORT& other) const throw (MEDEXCEPTION)
 828 {
 829   const char * LOC = "SUPPORT *MEDMEM::subtract(const SUPPORT * other) : ";
 830   BEGIN_OF(LOC);
 831   SUPPORT *ret;
 832   if (_entity!=other.getEntity())
 833     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Entities are different !"));
 834   if(!(*_mesh == *other.getMesh()))
 835     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Mesh are different !"));
 836   if(other._isOnAllElts)
 837     {
 838       ret=new SUPPORT;
 839       ret->setMesh(_mesh);
 840       ret->setEntity(_entity);
 841       return ret;
 842     }
 843   if(_isOnAllElts)
 844     return other.getComplement();
 845   int nbOfEltInThis=getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
 846   const int *nbsThis=_number->getValue();
 847   int nbOfEltInOther=other.getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
 848   const int *nbsOther=other._number->getValue();
 849   list<int> *ids=sub(nbsThis,nbOfEltInThis,nbsOther,nbOfEltInOther);
 850   if(_entity==MED_EN::MED_NODE)
 851     ret=_mesh->buildSupportOnNodeFromElementList(*ids,_entity);
 852   else
 853     ret=_mesh->buildSupportOnElementsFromElementList(*ids,_entity);
 854   delete ids;
 855   return ret;
 856   END_OF(LOC);
 857 }
 858
 859 /*!
 860   returns a new support the user has to delete. Entity is either MED_NODE to obtain node elements lying on boundary of "this"
 861   or MED_FACE,MED_EDGE (depends on the this->_mesh dimension).
 862  */
 863 SUPPORT *MEDMEM::SUPPORT::getBoundaryElements(MED_EN::medEntityMesh Entity) const throw (MEDEXCEPTION)
 864 {
 865   const char * LOC = "SUPPORT *MEDMEM::SUPPORT::getBoundaryElements(MED_EN::medEntityMesh Entity) : ";
 866   BEGIN_OF(LOC);
 867   int spaceDimension=_mesh->getSpaceDimension();
 868   MED_EN::medEntityMesh baseEntity=Entity;
 869   if (spaceDimension == 3)
 870     if (Entity!=MED_FACE)
 871       if(Entity==MED_NODE)
 872         baseEntity=MED_FACE;
 873       else
 874         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<Entity<<" !"));
 875   if (spaceDimension == 2)
 876     if (Entity!=MED_EDGE)
 877       if(Entity==MED_NODE)
 878         baseEntity=MED_EDGE;
 879       else
 880         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<Entity<<" !"));
 881   if(_isOnAllElts)
 882     return _mesh->getBoundaryElements(Entity);
 883
 884   const int * myConnectivityValue=_mesh->getReverseConnectivity(MED_DESCENDING);
 885   const int * myConnectivityIndex=_mesh->getReverseConnectivityIndex(MED_DESCENDING);
 886   int numberOf=_mesh->getNumberOfElements(baseEntity,MED_ALL_ELEMENTS);
 887   const int *ids=_number->getValue();
 888   set<int> idsSet(ids,ids+_totalNumberOfElements);
 889   list<int> myElementsList;
 890   for (int i=0;i<numberOf;i++)
 891     {
 892       int nbOfDP1EntitySharing=0;
 893       if(idsSet.find(myConnectivityValue[myConnectivityIndex[i]-1])!=idsSet.end())
 894         nbOfDP1EntitySharing++;
 895       if(idsSet.find(myConnectivityValue[myConnectivityIndex[i]])!=idsSet.end())
 896         nbOfDP1EntitySharing++;
 897       if(nbOfDP1EntitySharing==1)
 898         myElementsList.push_back(i+1);
 899     }
 900   if(Entity==MED_NODE)
 901     {
 902       return _mesh->buildSupportOnNodeFromElementList(myElementsList,baseEntity);
 903     }
 904   else
 905     {
 906       return _mesh->buildSupportOnElementsFromElementList(myElementsList,baseEntity);
 907     }
 908 }
 909
 910 /*!
 911   Method that fills this and updates all its attributes in order to lye on the the listOfNode.
 912  */
 913 void MEDMEM::SUPPORT::fillFromNodeList(const list<int>& listOfNode) throw (MEDEXCEPTION)
 914 {
 915   setEntity(MED_EN::MED_NODE);
 916   clearDataOnNumbers();
 917   int size=listOfNode.size();
 918   int totalNbInMesh=_mesh->getNumberOfElements(_entity,MED_ALL_ELEMENTS);
 919   if(totalNbInMesh==size)
 920     {
 921     _isOnAllElts=true;
 922     update();
 923     return;
 924     }
 925   else
 926     _isOnAllElts=false;
 927   int numberOfGeometricType=1;
 928   medGeometryElement* geometricType=new medGeometryElement[1];
 929   geometricType[0]=MED_NONE;
 930   int *numberOfGaussPoint=new int[1];
 931   numberOfGaussPoint[0]=1;
 932   int *numberOfElements=new int[1];
 933   numberOfElements[0]=size;
 934   int *mySkyLineArrayIndex=new int[2];
 935   mySkyLineArrayIndex[0]=1;
 936   mySkyLineArrayIndex[1]=1+numberOfElements[0];
 937   int *tab=new int[numberOfElements[0]];
 938   int i=0;
 939   for(list<int>::const_iterator iter2=listOfNode.begin();iter2!=listOfNode.end();iter2++)
 940     tab[i++]=*iter2;
 941   MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(1,numberOfElements[0],mySkyLineArrayIndex,tab,true);
 942   setNumberOfGeometricType(numberOfGeometricType);
 943   setGeometricType(geometricType);
 944   setNumberOfGaussPoint(numberOfGaussPoint);
 945   setNumberOfElements(numberOfElements);
 946   setTotalNumberOfElements(numberOfElements[0]);
 947   setNumber(mySkyLineArray);
 948
 949   delete[] numberOfElements;
 950   delete[] numberOfGaussPoint;
 951   delete[] geometricType;
 952 }
 953
 954 /*
 955   Method created to factorize code. This method fills the current SUPPORT on entity 'entity' containing all the entities contained in
 956   elements 'listOfElt' of entity 'entity'. Warning this method should be called after both the attributes this->_mesh and this->_entity are correctly set.
 957  */
 958 void MEDMEM::SUPPORT::fillFromElementList(const list<int>& listOfElt) throw (MEDEXCEPTION)
 959 {
 960   clearDataOnNumbers();
 961   int size=listOfElt.size();
 962   int totalNbInMesh=_mesh->getNumberOfElementsWithPoly(_entity,MED_ALL_ELEMENTS);
 963   if(totalNbInMesh==size){
 964     _isOnAllElts=true;
 965     update();
 966     return;
 967   }
 968   else
 969   _isOnAllElts=false;
 970   // Well, we must know how many geometric type we have found
 971   int * myListArray = new int[size] ;
 972   int id = 0 ;
 973   list<int>::const_iterator myElementsListIt ;
 974   for (myElementsListIt=listOfElt.begin();myElementsListIt!=listOfElt.end();myElementsListIt++)
 975     myListArray[id++]=(*myElementsListIt) ;
 976   int numberOfGeometricType ;
 977   medGeometryElement* geometricType ;
 978   int * numberOfGaussPoint ;
 979   int * numberOfElements ;
 980   int * mySkyLineArrayIndex ;
 981
 982   int numberOfType = _mesh->getNumberOfTypesWithPoly(_entity) ;
 983   if (numberOfType == 1) {
 984     numberOfGeometricType = 1 ;
 985     geometricType = new medGeometryElement[1] ;
 986     medGeometryElement *  allType = _mesh->getTypesWithPoly(_entity);
 987     geometricType[0] = allType[0] ;
 988     numberOfGaussPoint = new int[1] ;
 989     numberOfGaussPoint[0] = 1 ;
 990     numberOfElements = new int[1] ;
 991     numberOfElements[0] = size ;
 992     mySkyLineArrayIndex = new int[2] ;
 993     mySkyLineArrayIndex[0]=1 ;
 994     mySkyLineArrayIndex[1]=1+size ;
 995     delete [] allType;
 996   }
 997   else {// hemmm
 998     map<medGeometryElement,int> theType ;
 999     for (myElementsListIt=listOfElt.begin();myElementsListIt!=listOfElt.end();myElementsListIt++) {
1000       medGeometryElement myType = _mesh->getElementTypeWithPoly(_entity,*myElementsListIt) ;
1001       if (theType.find(myType) != theType.end() )
1002         theType[myType]+=1 ;
1003       else
1004         theType[myType]=1 ;
1005     }
1006     numberOfGeometricType = theType.size() ;
1007     geometricType = new medGeometryElement[numberOfGeometricType] ;
1008     numberOfGaussPoint = new int[numberOfGeometricType] ;
1009     numberOfElements = new int[numberOfGeometricType] ;
1010     mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
1011     int index = 0 ;
1012     mySkyLineArrayIndex[0]=1 ;
1013     map<medGeometryElement,int>::iterator theTypeIt ;
1014     for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
1015       geometricType[index] = (*theTypeIt).first ;
1016       numberOfGaussPoint[index] = 1 ;
1017       numberOfElements[index] = (*theTypeIt).second ;
1018       mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfElements[index] ;
1019       index++ ;
1020     }
1021   }
1022   MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray,true) ;
1023   setNumberOfGeometricType(numberOfGeometricType) ;
1024   setGeometricType(geometricType) ;
1025   setNumberOfGaussPoint(numberOfGaussPoint) ;
1026   setNumberOfElements(numberOfElements) ;
1027   setTotalNumberOfElements(size) ;
1028   setNumber(mySkyLineArray) ;
1029
1030   delete[] numberOfElements;
1031   delete[] numberOfGaussPoint;
1032   delete[] geometricType;
1033 }
1034
1035 /*! set the reference _mesh to Mesh */
1036 //--------------------------------------
1037 void SUPPORT::setMesh(MESH *Mesh)
1038 //--------------------------------------
1039 {
1040   if(_mesh)
1041     _mesh->removeReference();
1042   _mesh=Mesh;
1043   if(_mesh)
1044     _mesh->addReference();
1045 }
1046
1047 /*!
1048   addReference : reference counter presently disconnected in C++ -> just connected for client.
1049 */
1050 void MEDMEM::SUPPORT::addReference() const
1051 {
1052 }
1053
1054 /*!
1055   removeReference : reference counter presently disconnected in C++ -> just connected for client.
1056 */
1057 void MEDMEM::SUPPORT::removeReference() const
1058 {
1059 }