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 {
 442   const char * LOC="SUPPORT::intersecting(SUPPORT *) : ";
 443   BEGIN_OF(LOC);
 444   if (_entity!=mySupport->getEntity())
 445     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Entities are different !"));
 446   if(!(*_mesh == *mySupport->getMesh()))
 447     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Mesh are different !"));
 448   if(mySupport->isOnAllElements())
 449     return;
 450   if(_isOnAllElts)
 451     {
 452       *this=*mySupport;
 453       return;
 454     }
 455   if(_totalNumberOfElements==0)
 456     return;
 457   const int *ids=getNumber(MED_ALL_ELEMENTS);
 458   set<int> idsSet(ids,ids+getNumberOfElements(MED_ALL_ELEMENTS));
 459   const int *idsMySupport=mySupport->getNumber(MED_ALL_ELEMENTS);
 460   int mySupportSize=mySupport->getNumberOfElements(MED_ALL_ELEMENTS);
 461   set<int> idsSetMySupport(idsMySupport,idsMySupport+mySupportSize);
 462   set<int>::iterator iter;
 463   list<int> idsList;
 464   for(iter=idsSet.begin();iter!=idsSet.end();iter++)
 465     if(idsSetMySupport.find(*iter)!=idsSetMySupport.end())
 466       idsList.push_back(*iter);
 467   int size=idsSet.size();
 468   if(size!=0)
 469     {
 470       if(_entity==MED_NODE)
 471         fillFromNodeList(idsList);
 472       else
 473         fillFromElementList(idsList);
 474     }
 475   else
 476     {
 477       clearDataOnNumbers();
 478     }
 479   END_OF(LOC);
 480 };
 481
 482 /*!
 483   operator = perform et deep copy except for attribute _mesh
 484  */
 485 SUPPORT& MEDMEM::SUPPORT::operator=(const SUPPORT &other)
 486 {
 487   const char * LOC = "SUPPORT::operator=(const SUPPORT & m) : ";
 488   BEGIN_OF(LOC);
 489
 490   _name=other._name;
 491   _description=other._description;
 492   _mesh=other._mesh;
 493   _entity=other._entity;
 494   _isOnAllElts=other._isOnAllElts;
 495   clearDataOnNumbers();
 496   _numberOfGeometricType=other._numberOfGeometricType;
 497   _totalNumberOfElements=other._totalNumberOfElements;
 498
 499   if(other._geometricType != NULL)
 500     {
 501       _geometricType=new medGeometryElement[other._numberOfGeometricType];
 502       memcpy(_geometricType,other._geometricType,other._numberOfGeometricType*sizeof(medGeometryElement));
 503     }
 504   if (other._numberOfGaussPoint != NULL)
 505     {
 506       _numberOfGaussPoint=new int[other._numberOfGeometricType];
 507       memcpy(_numberOfGaussPoint,other._numberOfGaussPoint,other._numberOfGeometricType*sizeof(int));
 508     }
 509   if(other._numberOfElements != NULL)
 510     {
 511       _numberOfElements=new int[_numberOfGeometricType];
 512       memcpy(_numberOfElements,other._numberOfElements,_numberOfGeometricType*sizeof(int));
 513     }
 514   if(!_isOnAllElts)
 515     {
 516       _number=new MEDSKYLINEARRAY(* other._number);
 517     }
 518   return *this;
 519 }
 520
 521 /*!
 522   Method that cleans up all the fields related to _numbers. Defined for code factorization.
 523  */
 524 //--------------------------------------------------
 525 void MEDMEM::SUPPORT::clearDataOnNumbers()
 526 //--------------------------------------------------
 527 {
 528   _numberOfGeometricType=0;
 529   _totalNumberOfElements=0;
 530   if(_geometricType)
 531     {
 532       delete [] _geometricType;
 533       _geometricType=(medGeometryElement *) NULL;
 534     }
 535   if(_numberOfGaussPoint)
 536     {
 537       delete [] _numberOfGaussPoint;
 538       _numberOfGaussPoint=(int *) NULL;
 539     }
 540   if(_numberOfElements)
 541     {
 542       delete [] _numberOfElements;
 543       _numberOfElements=(int *) NULL;
 544     }
 545   if(_number)
 546     {
 547       delete _number;
 548       _number=(MEDSKYLINEARRAY *) NULL;
 549     }
 550 }
 551
 552 /*!
 553   operator == This operator does not compare attributs _name and _description.
 554 */
 555 //--------------------------------------------------
 556 bool MEDMEM::SUPPORT::operator == (const SUPPORT &support) const
 557 //--------------------------------------------------
 558 {
 559   const char * LOC = "bool SUPPORT::operator ==(const SUPPORT &support) const : ";
 560
 561   BEGIN_OF(LOC);
 562
 563   bool operatorReturn = false;
 564
 565   operatorReturn = (*_mesh == *support._mesh) && (_entity == support._entity) &&
 566     (_numberOfGeometricType == support._numberOfGeometricType) &&
 567     ((_isOnAllElts && support._isOnAllElts) || (!_isOnAllElts && !support._isOnAllElts)) &&
 568     (_totalNumberOfElements == support._totalNumberOfElements);
 569
 570   if (operatorReturn)
 571     {
 572       if (!_isOnAllElts)
 573         {
 574           for (int i=0; i<_numberOfGeometricType; i++)
 575             {
 576               operatorReturn = operatorReturn &&
 577                 (_geometricType[i] == support._geometricType[i]) &&
 578                 (_numberOfElements[i] == support._numberOfElements[i]) &&
 579                 (_numberOfGaussPoint[i] == support._numberOfGaussPoint[i]);
 580
 581               if (operatorReturn)
 582                 {
 583                   for (int j=0; j<_numberOfElements[i]; j++)
 584                     {
 585                       operatorReturn = operatorReturn &&
 586                         (getNumber(_geometricType[i])[j] ==
 587                          support.getNumber(_geometricType[i])[j]);
 588                     }
 589                 }
 590             }
 591         }
 592     }
 593
 594   END_OF(LOC);
 595
 596   return operatorReturn;
 597 };
 598
 599 void SUPPORT::changeElementsNbs(MED_EN::medEntityMesh entity, const int *renumberingFromOldToNew, int limitNbClassicPoly, const int *renumberingFromOldToNewPoly)
 600 {
 601   if(entity != _entity)
 602     throw MEDEXCEPTION("SUPPORT::changeElementsNbs : Renumbering on a mismatch entity");
 603   list<int> newNbs;
 604   if(!_isOnAllElts)
 605     {
 606       const int *oldNbs=_number->getValue();
 607       for(int i=0;i<_totalNumberOfElements;i++)
 608         {
 609           int globNb=oldNbs[i];
 610           if(globNb<=limitNbClassicPoly)
 611             newNbs.push_back(renumberingFromOldToNew[globNb-1]);
 612           else
 613             newNbs.push_back(renumberingFromOldToNewPoly[globNb-limitNbClassicPoly-1]);
 614         }
 615       newNbs.sort();
 616       fillFromElementList(newNbs);
 617     }
 618   else
 619     update();
 620 }
 621
 622 /*!
 623   operator == + in case false a test if coordinates and connectivity of _mesh and support->_mesh are the same
 624 */
 625 bool MEDMEM::SUPPORT::deepCompare(const SUPPORT &support) const
 626 {
 627   bool operatorReturn =(_entity == support._entity) &&
 628     (_numberOfGeometricType == support._numberOfGeometricType) &&
 629     ( (_isOnAllElts && support._isOnAllElts) || (!_isOnAllElts  && !support._isOnAllElts) ) &&
 630     (_totalNumberOfElements == support._totalNumberOfElements);
 631   if (operatorReturn)
 632     {
 633       if (!_isOnAllElts)
 634         {
 635           for (int i=0; i<_numberOfGeometricType && operatorReturn; i++)
 636             {
 637               operatorReturn = (_geometricType[i] == support._geometricType[i]) &&
 638                 (_numberOfElements[i] == support._numberOfElements[i]) &&
 639                 (_numberOfGaussPoint[i] == support._numberOfGaussPoint[i]);
 640               if (operatorReturn)
 641                 {
 642                   for (int j=0; j<_numberOfElements[i]; j++)
 643                     {
 644                       operatorReturn = (getNumber(_geometricType[i])[j] ==
 645                          support.getNumber(_geometricType[i])[j]);
 646                     }
 647                 }
 648             }
 649         }
 650     }
 651   if(operatorReturn)
 652     {
 653       if(!(*_mesh == *support._mesh))
 654         {
 655           return _mesh->deepCompare(*support._mesh);
 656         }
 657     }
 658   return operatorReturn;
 659 }
 660
 661 /*!
 662  States if this is included in other.
 663  */
 664 bool MEDMEM::SUPPORT::belongsTo(const SUPPORT& other, bool deepCompare) const
 665 {
 666   if(!(*_mesh == *other._mesh))
 667     {
 668       if(!deepCompare)
 669         return false;
 670       if(!_mesh->deepCompare(*other._mesh))
 671         return false;
 672     }
 673   if(_entity!=other._entity)
 674     return false;
 675   if(other._isOnAllElts)
 676     return true;
 677   if(_isOnAllElts && !other._isOnAllElts)
 678     return false;
 679   if(_numberOfGeometricType>other._numberOfGeometricType)
 680     return false;
 681   for(int i=0; i<_numberOfGeometricType; i++)
 682     {
 683       MED_EN::medGeometryElement curGeomType=_geometricType[i];
 684       int iOther=-1;
 685       for(int j=0; j<other._numberOfGeometricType; j++)
 686           if(other._geometricType[j]==curGeomType)
 687               iOther=j;
 688       if(iOther==-1)
 689         return false;
 690       if(_numberOfElements[i]>other._numberOfElements[iOther])
 691         return false;
 692       const int *numbers1=_number->getI(i+1);
 693       const int *numbers2=other._number->getI(iOther+1);
 694       for (int k=0; k<_numberOfElements[i]; k++)
 695         {
 696           bool found=false;
 697           for(int l=0;l<other._numberOfElements[iOther] && !found;l++)
 698             {
 699               if(numbers1[k]==numbers2[l])
 700                 found=true;
 701             }
 702           if(!found)
 703             return false;
 704         }
 705     }
 706   return true;
 707 }
 708 /*!
 709   Method used to sort array of id.
 710  */
 711 int compareId(const void *x, const void *y)
 712 {
 713   const int *x1=(const int *)x;
 714   const int *y1=(const int *)y;
 715   if(*x1<*y1)
 716     return -1;
 717   else if(*x1>*y1)
 718     return 1;
 719   else
 720     return 0;
 721 }
 722
 723 /*!
 724   performs a common operation : Sub builds a sorted int array that is obtained by supression of all ids contained
 725   in array defined by (idsToSuppress,lgthIdsToSuppress) from array [start ... end]
 726   Example sub(0,7,{1,2,5},3) => {0,3,4,6,7} - WARNING returned list should be deallocated !
 727  */
 728 list<int> *MEDMEM::SUPPORT::sub(int start,int end,const int *idsToSuppress,int lgthIdsToSuppress)
 729 {
 730   int size=end-start+1;
 731   int sizeRet=size-lgthIdsToSuppress;
 732   list<int> *ret;
 733   if(sizeRet<0)
 734     throw MEDEXCEPTION("MEDMEM::SUPPORT::sub");
 735   else if(sizeRet==0)
 736     {
 737       return 0;
 738     }
 739   if(idsToSuppress==0)
 740     {
 741       ret=new list<int>;
 742       for(int l=0;l<size;l++)
 743         ret->push_back(start+l);
 744       return ret;
 745     }
 746   ret=new list<int>;
 747   int *temp=new int[lgthIdsToSuppress];
 748   memcpy(temp,idsToSuppress,sizeof(int)*lgthIdsToSuppress);
 749   qsort(temp,lgthIdsToSuppress,sizeof(int),compareId);
 750   int k=0;
 751   for(int i=start;i<=end;i++)
 752     if(temp[k]!=i)
 753       ret->push_back(i);
 754     else
 755       k++;
 756   delete [] temp;
 757   return ret;
 758 }
 759
 760 /*!
 761   performs a common operation : Sub builds a sorted int array that is obtained by supression of all ids contained
 762   in array defined by (idsToSuppress,lgthIdsToSuppress) from array [start ... end]
 763   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 !
 764  */
 765 list<int> *MEDMEM::SUPPORT::sub(const int *ids,int lgthIds,const int *idsToSuppress,int lgthIdsToSuppress)
 766 {
 767   list<int> *ret;
 768   int i,j=0;
 769   if(lgthIds<0)
 770     throw MEDEXCEPTION("MEDMEM::SUPPORT::sub");
 771   else if(lgthIds==0)
 772     return 0;
 773   ret=new list<int>;
 774   int *temp1=new int[lgthIds];
 775   memcpy(temp1,ids,sizeof(int)*lgthIds);
 776   qsort(temp1,lgthIds,sizeof(int),compareId);
 777   int *temp2=new int[lgthIdsToSuppress];
 778   memcpy(temp2,idsToSuppress,sizeof(int)*lgthIdsToSuppress);
 779   qsort(temp2,lgthIdsToSuppress,sizeof(int),compareId);
 780   for(i=0;i<lgthIds;)
 781     {
 782       if(j>=lgthIdsToSuppress)
 783           ret->push_back(temp1[i++]);
 784       else if(temp1[i]>temp2[j])
 785         j++;
 786       else if(temp1[i]<temp2[j])
 787         ret->push_back(temp1[i++]);
 788       else
 789         i++;
 790     }
 791   delete [] temp1;
 792   delete [] temp2;
 793   return ret;
 794 }
 795
 796 /*!
 797   returns a new SUPPORT (responsability to caller to destroy it)
 798   that is the complement to "this", lying on the same entity than "this".
 799  */
 800 SUPPORT *MEDMEM::SUPPORT::getComplement() const
 801 {
 802   SUPPORT *ret;
 803   const int nbOfElt=_mesh->getNumberOfElements(_entity,MED_EN::MED_ALL_ELEMENTS);
 804   int nbOfEltInSupp=getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
 805   if(_isOnAllElts || nbOfElt==nbOfEltInSupp)
 806     {
 807       ret=new SUPPORT;
 808       ret->setMesh(_mesh);
 809       ret->setEntity(_entity);
 810       string name="Complement of ";
 811       name+=_name;
 812       ret->setName(name);
 813       return ret;
 814     }
 815   const int *nbs=_number->getValue();
 816   list<int> *ids=sub(1,nbOfElt,nbs,nbOfEltInSupp);
 817   if(_entity==MED_EN::MED_NODE)
 818     ret=_mesh->buildSupportOnNodeFromElementList(*ids,_entity);
 819   else
 820     ret=_mesh->buildSupportOnElementsFromElementList(*ids,_entity);
 821   delete ids;
 822   return ret;
 823 }
 824
 825 /*!
 826   returns a new support the user should delete.
 827  */
 828 SUPPORT *MEDMEM::SUPPORT::substract(const SUPPORT& other) const throw (MEDEXCEPTION)
 829 {
 830   const char * LOC = "SUPPORT *MEDMEM::subtract(const SUPPORT * other) : ";
 831   BEGIN_OF(LOC);
 832   SUPPORT *ret;
 833   if (_entity!=other.getEntity())
 834     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Entities are different !"));
 835   if(!(*_mesh == *other.getMesh()))
 836     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Mesh are different !"));
 837   if(other._isOnAllElts)
 838     {
 839       ret=new SUPPORT;
 840       ret->setMesh(_mesh);
 841       ret->setEntity(_entity);
 842       return ret;
 843     }
 844   if(_isOnAllElts)
 845     return other.getComplement();
 846   int nbOfEltInThis=getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
 847   const int *nbsThis=_number->getValue();
 848   int nbOfEltInOther=other.getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
 849   const int *nbsOther=other._number->getValue();
 850   list<int> *ids=sub(nbsThis,nbOfEltInThis,nbsOther,nbOfEltInOther);
 851   if(_entity==MED_EN::MED_NODE)
 852     ret=_mesh->buildSupportOnNodeFromElementList(*ids,_entity);
 853   else
 854     ret=_mesh->buildSupportOnElementsFromElementList(*ids,_entity);
 855   delete ids;
 856   return ret;
 857   END_OF(LOC);
 858 }
 859
 860 /*!
 861   returns a new support the user has to delete. Entity is either MED_NODE to obtain node elements lying on boundary of "this"
 862   or MED_FACE,MED_EDGE (depends on the this->_mesh dimension).
 863  */
 864 SUPPORT *MEDMEM::SUPPORT::getBoundaryElements(MED_EN::medEntityMesh Entity) const throw (MEDEXCEPTION)
 865 {
 866   const char * LOC = "SUPPORT *MEDMEM::SUPPORT::getBoundaryElements(MED_EN::medEntityMesh Entity) : ";
 867   BEGIN_OF(LOC);
 868   int spaceDimension=_mesh->getSpaceDimension();
 869   MED_EN::medEntityMesh baseEntity=Entity;
 870   if (spaceDimension == 3)
 871     if (Entity!=MED_FACE)
 872       if(Entity==MED_NODE)
 873         baseEntity=MED_FACE;
 874       else
 875         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<Entity<<" !"));
 876   if (spaceDimension == 2)
 877     if (Entity!=MED_EDGE)
 878       if(Entity==MED_NODE)
 879         baseEntity=MED_EDGE;
 880       else
 881         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<Entity<<" !"));
 882   if(_isOnAllElts)
 883     return _mesh->getBoundaryElements(Entity);
 884
 885   const int * myConnectivityValue=_mesh->getReverseConnectivity(MED_DESCENDING);
 886   const int * myConnectivityIndex=_mesh->getReverseConnectivityIndex(MED_DESCENDING);
 887   int numberOf=_mesh->getNumberOfElements(baseEntity,MED_ALL_ELEMENTS);
 888   const int *ids=_number->getValue();
 889   set<int> idsSet(ids,ids+_totalNumberOfElements);
 890   list<int> myElementsList;
 891   for (int i=0;i<numberOf;i++)
 892     {
 893       int nbOfDP1EntitySharing=0;
 894       if(idsSet.find(myConnectivityValue[myConnectivityIndex[i]-1])!=idsSet.end())
 895         nbOfDP1EntitySharing++;
 896       if(idsSet.find(myConnectivityValue[myConnectivityIndex[i]])!=idsSet.end())
 897         nbOfDP1EntitySharing++;
 898       if(nbOfDP1EntitySharing==1)
 899         myElementsList.push_back(i+1);
 900     }
 901   if(Entity==MED_NODE)
 902     {
 903       return _mesh->buildSupportOnNodeFromElementList(myElementsList,baseEntity);
 904     }
 905   else
 906     {
 907       return _mesh->buildSupportOnElementsFromElementList(myElementsList,baseEntity);
 908     }
 909 }
 910
 911 /*!
 912   Method that fills this and updates all its attributes in order to lye on the the listOfNode.
 913  */
 914 void MEDMEM::SUPPORT::fillFromNodeList(const list<int>& listOfNode) throw (MEDEXCEPTION)
 915 {
 916   setEntity(MED_EN::MED_NODE);
 917   clearDataOnNumbers();
 918   int size=listOfNode.size();
 919   int totalNbInMesh=_mesh->getNumberOfElements(_entity,MED_ALL_ELEMENTS);
 920   if(totalNbInMesh==size)
 921     {
 922     _isOnAllElts=true;
 923     update();
 924     return;
 925     }
 926   else
 927     _isOnAllElts=false;
 928   int numberOfGeometricType=1;
 929   medGeometryElement* geometricType=new medGeometryElement[1];
 930   geometricType[0]=MED_NONE;
 931   int *numberOfGaussPoint=new int[1];
 932   numberOfGaussPoint[0]=1;
 933   int *numberOfElements=new int[1];
 934   numberOfElements[0]=size;
 935   int *mySkyLineArrayIndex=new int[2];
 936   mySkyLineArrayIndex[0]=1;
 937   mySkyLineArrayIndex[1]=1+numberOfElements[0];
 938   int *tab=new int[numberOfElements[0]];
 939   int i=0;
 940   for(list<int>::const_iterator iter2=listOfNode.begin();iter2!=listOfNode.end();iter2++)
 941     tab[i++]=*iter2;
 942   MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(1,numberOfElements[0],mySkyLineArrayIndex,tab,true);
 943   setNumberOfGeometricType(numberOfGeometricType);
 944   setGeometricType(geometricType);
 945   setNumberOfGaussPoint(numberOfGaussPoint);
 946   setNumberOfElements(numberOfElements);
 947   setTotalNumberOfElements(numberOfElements[0]);
 948   setNumber(mySkyLineArray);
 949
 950   delete[] numberOfElements;
 951   delete[] numberOfGaussPoint;
 952   delete[] geometricType;
 953 }
 954
 955 /*
 956   Method created to factorize code. This method fills the current SUPPORT on entity 'entity' containing all the entities contained in
 957   elements 'listOfElt' of entity 'entity'. Warning this method should be called after both the attributes this->_mesh and this->_entity are correctly set.
 958  */
 959 void MEDMEM::SUPPORT::fillFromElementList(const list<int>& listOfElt) throw (MEDEXCEPTION)
 960 {
 961   clearDataOnNumbers();
 962   int size=listOfElt.size();
 963   int totalNbInMesh=_mesh->getNumberOfElementsWithPoly(_entity,MED_ALL_ELEMENTS);
 964   if(totalNbInMesh==size){
 965     _isOnAllElts=true;
 966     update();
 967     return;
 968   }
 969   else
 970   _isOnAllElts=false;
 971   // Well, we must know how many geometric type we have found
 972   int * myListArray = new int[size] ;
 973   int id = 0 ;
 974   list<int>::const_iterator myElementsListIt ;
 975   for (myElementsListIt=listOfElt.begin();myElementsListIt!=listOfElt.end();myElementsListIt++)
 976     myListArray[id++]=(*myElementsListIt) ;
 977   int numberOfGeometricType ;
 978   medGeometryElement* geometricType ;
 979   int * numberOfGaussPoint ;
 980   int * numberOfElements ;
 981   int * mySkyLineArrayIndex ;
 982
 983   int numberOfType = _mesh->getNumberOfTypesWithPoly(_entity) ;
 984   if (numberOfType == 1) {
 985     numberOfGeometricType = 1 ;
 986     geometricType = new medGeometryElement[1] ;
 987     medGeometryElement *  allType = _mesh->getTypesWithPoly(_entity);
 988     geometricType[0] = allType[0] ;
 989     numberOfGaussPoint = new int[1] ;
 990     numberOfGaussPoint[0] = 1 ;
 991     numberOfElements = new int[1] ;
 992     numberOfElements[0] = size ;
 993     mySkyLineArrayIndex = new int[2] ;
 994     mySkyLineArrayIndex[0]=1 ;
 995     mySkyLineArrayIndex[1]=1+size ;
 996     delete [] allType;
 997   }
 998   else {// hemmm
 999     map<medGeometryElement,int> theType ;
1000     for (myElementsListIt=listOfElt.begin();myElementsListIt!=listOfElt.end();myElementsListIt++) {
1001       medGeometryElement myType = _mesh->getElementTypeWithPoly(_entity,*myElementsListIt) ;
1002       if (theType.find(myType) != theType.end() )
1003         theType[myType]+=1 ;
1004       else
1005         theType[myType]=1 ;
1006     }
1007     numberOfGeometricType = theType.size() ;
1008     geometricType = new medGeometryElement[numberOfGeometricType] ;
1009     numberOfGaussPoint = new int[numberOfGeometricType] ;
1010     numberOfElements = new int[numberOfGeometricType] ;
1011     mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
1012     int index = 0 ;
1013     mySkyLineArrayIndex[0]=1 ;
1014     map<medGeometryElement,int>::iterator theTypeIt ;
1015     for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
1016       geometricType[index] = (*theTypeIt).first ;
1017       numberOfGaussPoint[index] = 1 ;
1018       numberOfElements[index] = (*theTypeIt).second ;
1019       mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfElements[index] ;
1020       index++ ;
1021     }
1022   }
1023   MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray,true) ;
1024   setNumberOfGeometricType(numberOfGeometricType) ;
1025   setGeometricType(geometricType) ;
1026   setNumberOfGaussPoint(numberOfGaussPoint) ;
1027   setNumberOfElements(numberOfElements) ;
1028   setTotalNumberOfElements(size) ;
1029   setNumber(mySkyLineArray) ;
1030
1031   delete[] numberOfElements;
1032   delete[] numberOfGaussPoint;
1033   delete[] geometricType;
1034 }
1035
1036 /*! set the reference _mesh to Mesh */
1037 //--------------------------------------
1038 void SUPPORT::setMesh(MESH *Mesh)
1039 //--------------------------------------
1040 {
1041   if(_mesh)
1042     _mesh->removeReference();
1043   _mesh=Mesh;
1044   if(_mesh)
1045     _mesh->addReference();
1046 }
1047
1048 /*!
1049   addReference : reference counter presently disconnected in C++ -> just connected for client.
1050 */
1051 void MEDMEM::SUPPORT::addReference() const
1052 {
1053 }
1054
1055 /*!
1056   removeReference : reference counter presently disconnected in C++ -> just connected for client.
1057 */
1058 void MEDMEM::SUPPORT::removeReference() const
1059 {
1060 }