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