src/MEDMEM/MEDMEM_Support.cxx

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