src/MEDMEM/MEDMEM_Support.cxx

   1 using namespace std;
   2 /*
   3  File Support.cxx
   4  $Header$
   5 */
   6
   7 #include <set>
   8 #include <algorithm>
   9 #include <list>
  10
  11 //#include "utilities.h"
  12 //#include "MEDMEM_define.hxx"
  13 #include "MEDMEM_DriversDef.hxx"
  14 #include "MEDMEM_Support.hxx"
  15 //#include "MEDMEM_Family.hxx"
  16 //#include "MEDMEM_Group.hxx"
  17 #include "MEDMEM_Mesh.hxx"
  18
  19 using namespace MED_EN;
  20 using namespace MEDMEM;
  21
  22 /* This class is a generic class for family and group */
  23
  24 /*!
  25   Constructor.
  26 */
  27 //--------------------------------------------------------------------------
  28 SUPPORT::SUPPORT(): _name(""),  _description(""), _mesh((MESH*)NULL),
  29                     _entity(MED_CELL), _numberOfGeometricType(0),
  30                     _geometricType((medGeometryElement*)NULL),
  31                     _numberOfGaussPoint((int*)NULL),
  32                     //_geometricTypeNumber((int*)NULL),
  33                     _isOnAllElts(false),
  34                     _numberOfElements((int*)NULL),
  35                     _totalNumberOfElements(0),
  36                     _number((MEDSKYLINEARRAY*)NULL)
  37 //--------------------------------------------------------------------------
  38 {
  39     MESSAGE("SUPPORT::SUPPORT()");
  40 };
  41
  42 /*!
  43   Constructor.
  44 */
  45 //--------------------------------------------------------------------------
  46 SUPPORT::SUPPORT(MESH* Mesh, string Name/*=""*/, medEntityMesh Entity/*=MED_CELL*/):
  47                 _name(Name), _description(""), _mesh(Mesh), _entity(Entity),
  48                 _numberOfGeometricType(0),
  49                 _geometricType((medGeometryElement*)NULL),
  50                 _numberOfGaussPoint((int*)NULL),
  51                 //_geometricTypeNumber((int*)NULL),
  52                 _isOnAllElts(true),
  53                 _numberOfElements((int*)NULL),
  54                 _totalNumberOfElements(0),
  55                 _number((MEDSKYLINEARRAY*)NULL)
  56 //--------------------------------------------------------------------------
  57 {
  58   MESSAGE("SUPPORT::SUPPORT(MESH*Mesh,string Name,medEntityMesh Entity)");
  59   update() ;
  60 };
  61
  62 /*!
  63   Copy constructor.
  64 */
  65 //--------------------------------------------------------------------------
  66 SUPPORT::SUPPORT(const SUPPORT & m)
  67 //--------------------------------------------------------------------------
  68 {
  69   const char * LOC = "SUPPORT::SUPPORT(SUPPORT & m) : " ;
  70   BEGIN_OF(LOC) ;
  71
  72   _name = m._name;
  73   _description = m._description;
  74   _mesh = m._mesh; // on recopie uniquement l'adresse
  75   _entity = m._entity;
  76   _numberOfGeometricType = m._numberOfGeometricType;
  77   if (m._geometricType != NULL)
  78     {
  79       _geometricType = new medGeometryElement[m._numberOfGeometricType];
  80       memcpy(_geometricType,m._geometricType,m._numberOfGeometricType*sizeof(medGeometryElement));
  81     }
  82   else
  83     _geometricType = (medGeometryElement *) NULL;
  84   if (m._numberOfGaussPoint != NULL)
  85     {
  86       _numberOfGaussPoint = new int[m._numberOfGeometricType];
  87       memcpy(_numberOfGaussPoint,m._numberOfGaussPoint,m._numberOfGeometricType*sizeof(int));
  88     }
  89   else
  90     _numberOfGaussPoint = (int *) NULL;
  91 //    if (m._geometricTypeNumber != NULL)
  92 //      {
  93 //        _geometricTypeNumber = new int[m._numberOfGeometricType];
  94 //        memcpy(_geometricTypeNumber,m._geometricTypeNumber,m._numberOfGeometricType*sizeof(int));
  95 //      }
  96 //    else
  97 //      _geometricTypeNumber = (int *) NULL;
  98   _isOnAllElts = m._isOnAllElts;
  99   if (m._numberOfElements != NULL)
 100     {
 101       _numberOfElements = new int[_numberOfGeometricType];
 102       memcpy(_numberOfElements,m._numberOfElements,_numberOfGeometricType*sizeof(int));
 103     }
 104   else
 105     _numberOfElements = (int *) NULL;
 106   _totalNumberOfElements = m._totalNumberOfElements;
 107   if (m._isOnAllElts == false)
 108     _number = new MEDSKYLINEARRAY(* m._number);
 109   else
 110     _number = (MEDSKYLINEARRAY *) NULL;
 111
 112   END_OF(LOC) ;
 113 };
 114
 115
 116 /*!
 117   Destructor.
 118 */
 119 //-----------------
 120 SUPPORT::~SUPPORT()
 121 //-----------------
 122 {
 123   MESSAGE("Destructeur ~SUPPORT()");
 124   if (_geometricType != (medGeometryElement *) NULL)
 125     delete [] _geometricType ;
 126   if (_numberOfGaussPoint != (int *) NULL)
 127     delete [] _numberOfGaussPoint ;
 128   //      if (_geometricTypeNumber!=NULL)
 129   //                    delete[] _geometricTypeNumber ;
 130   if (_numberOfElements != (int *) NULL)
 131     delete[] _numberOfElements ;
 132   if (_number != (MEDSKYLINEARRAY *) NULL)
 133     delete _number ;
 134 }
 135
 136 /*!
 137   operator <<.
 138 */
 139 //--------------------------------------------------
 140 ostream & MEDMEM::operator<<(ostream &os, const SUPPORT &my)
 141 //--------------------------------------------------
 142 {
 143   os << "Name : "<< my.getName() << endl ;
 144   os << "Description : "<< my.getDescription() << endl ;
 145   os << "Mesh name : ";
 146   if (my.getMesh() == NULL)
 147     os << "ERROR : Mesh not defined !" << endl ;
 148   else
 149     os << my._mesh->getName() << endl ;
 150   os << "Entity : "<< my._entity << endl;
 151   os << "Entity list : "<< endl;
 152   if (!(my._isOnAllElts)) {
 153     int numberoftypes = my._numberOfGeometricType ;
 154     os << "NumberOfTypes : "<<numberoftypes<<endl;
 155     medGeometryElement * types = my._geometricType;
 156     for (int j=0;j<numberoftypes;j++) {
 157       int numberOfElements = my._numberOfElements[j];
 158       os << "    * Type "<<types[j]<<" : there is(are) "<<numberOfElements<<" element(s) : ";
 159       const int * number = my.getNumber(types[j]);
 160       SCRUTE(number);
 161       for (int k=0; k<numberOfElements;k++)
 162         os << number[k] << " ";
 163       os << endl ;
 164     }
 165   } else
 166     os << "Is on all entities !"<< endl;
 167
 168   return os ;
 169 }
 170
 171 /*!
 172   Updade the SUPPORT attributs with rigth MESH information.
 173
 174   It has an effect only if SUPPORT is on all elements.
 175
 176   No more need in future release.
 177 */
 178 //-------------------
 179 void SUPPORT::update()
 180 //-------------------
 181 {
 182   const char * LOC = "SUPPORT::update() : " ;
 183   BEGIN_OF(LOC) ;
 184
 185   if (_isOnAllElts) {
 186     if (_entity == MED_NODE) {
 187       _numberOfGeometricType=1 ;
 188       _geometricType=new medGeometryElement[1] ;
 189       _geometricType[0]=MED_NONE ;
 190       _numberOfElements = new int[1] ;
 191       _numberOfElements[0]=_mesh->getNumberOfNodes();
 192       _totalNumberOfElements=_numberOfElements[0];
 193       _numberOfGaussPoint = new int[1] ;
 194       _numberOfGaussPoint[0]=1;
 195     } else { // we duplicate information from _mesh
 196       _numberOfGeometricType=_mesh->getNumberOfTypes(_entity);
 197       if (_geometricType == (medGeometryElement *) NULL)
 198         _geometricType = new medGeometryElement[_numberOfGeometricType] ;
 199       memcpy(_geometricType,_mesh->getTypes(_entity),_numberOfGeometricType*sizeof(medGeometryElement));
 200       if (_numberOfElements == (int *) NULL)
 201         _numberOfElements = new int[_numberOfGeometricType] ;
 202       if (_numberOfGaussPoint == (int *) NULL)
 203         _numberOfGaussPoint = new int[_numberOfGeometricType] ;
 204       _totalNumberOfElements=0;
 205       for (int i=0;i<_numberOfGeometricType;i++) {
 206         _numberOfElements[i]=_mesh->getNumberOfElements(_entity,_geometricType[i]) ;
 207         _totalNumberOfElements+=_numberOfElements[i];
 208         _numberOfGaussPoint[i]=1 ;
 209       }
 210     }
 211   }
 212   END_OF(LOC);
 213 };
 214
 215 /*!
 216   Blend the given SUPPORT mySupport into the calling object SUPPORT.
 217 */
 218 //-------------------
 219 void SUPPORT::blending(SUPPORT * mySupport) throw (MEDEXCEPTION)
 220 //-------------------
 221 {
 222   const char * LOC = "SUPPORT::blending() : " ;
 223   BEGIN_OF(LOC) ;
 224
 225   MESSAGE(LOC<< "SUPPORT entry : " << *mySupport) ;
 226
 227   // on same entity :
 228   if ( _entity != mySupport->getEntity() )
 229     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Entities are different !"));
 230
 231   int * tmp_NumberOfElementsInType = new int[MED_NBR_GEOMETRIE_MAILLE];
 232   medGeometryElement * myType = new medGeometryElement[MED_NBR_GEOMETRIE_MAILLE];
 233   int * whereIsType = new int[MED_NBR_GEOMETRIE_MAILLE];
 234   //MESH_ENTITIES myMeshEntities() ;
 235   list<MED_FR::med_geometrie_element>::const_iterator listIt ;
 236   int it=0 ;
 237   for(listIt=(MED_FR::meshEntities[(MED_FR::med_entite_maillage)_entity]).begin();listIt!=(MED_FR::meshEntities[(MED_FR::med_entite_maillage)_entity]).end();listIt++) {
 238     tmp_NumberOfElementsInType[it]=0;
 239     whereIsType[it]=0 ;
 240     try {
 241       int tmp_int = 0;
 242       tmp_int = getNumberOfElements((medGeometryElement)(*listIt)) ;
 243       tmp_NumberOfElementsInType[it]+=tmp_int ;
 244       whereIsType[it]+=1 ;
 245     }
 246     catch (MEDEXCEPTION & ex) { SCRUTE(sizeof(ex)); };
 247     try {
 248       int tmp_int = 0;
 249       tmp_int = mySupport->getNumberOfElements((medGeometryElement)(*listIt)) ;
 250       tmp_NumberOfElementsInType[it]+=tmp_int ;
 251       whereIsType[it]+=2 ;
 252     }
 253     catch (const MEDEXCEPTION & ex) {};
 254     if (whereIsType[it]!=0) {
 255       myType[it]=(medGeometryElement)(*listIt) ;
 256       SCRUTE(myType[it]);SCRUTE(it);SCRUTE((*listIt));
 257       it++;
 258     }
 259   }
 260   // set new value :
 261 //   int * numberOfElements=_numberOfElements ;
 262 //   _numberOfElements = new int[it] ;
 263   int * numberOfElements= new int[it];
 264   _totalNumberOfElements = 0 ;
 265   //int totalSize = 0 ;
 266   int ** tmp_array = new (int*)[it];
 267   for (int i=0;i<it;i++) {
 268     int numberOfElementsInType = tmp_NumberOfElementsInType[i] ;
 269     numberOfElements[i] = numberOfElementsInType ;
 270     tmp_array[i] = new int[numberOfElementsInType] ;
 271     //totalSize+=numberOfElementsInType*(myType[i]%100) ;
 272     _totalNumberOfElements+=numberOfElementsInType ;
 273     if (whereIsType[i] == 1) { // only first Support
 274       memcpy(tmp_array[i],getNumber(myType[i]),sizeof(int)*numberOfElementsInType);
 275     } else if (whereIsType[i] == 2) { // only second Support
 276       memcpy(tmp_array[i],mySupport->getNumber(myType[i]),sizeof(int)*numberOfElementsInType);
 277     } else if (whereIsType[i] == 3) { // more difficult :-)
 278       set<int> elementList ;
 279       //int i1 = 0 ; !! UNUSED VARIABLE !!
 280       //int i2 = 0 ; !!UNUSED VARIABLE !!
 281       int ii = 0 ;
 282       const int * number1 = getNumber(myType[i]) ;
 283       const int * number2 = mySupport->getNumber(myType[i]) ;
 284
 285       SCRUTE(number1);
 286       SCRUTE(number2);
 287
 288       int numberOfElements1 = getNumberOfElements(myType[i]) ;
 289       int numberOfElements2 = mySupport->getNumberOfElements(myType[i]) ;
 290
 291       SCRUTE(numberOfElements1);
 292       SCRUTE(numberOfElements2);
 293
 294       MESSAGE(LOC << " Type : " << myType[i] << " " << i);
 295
 296       for(int j=0;j<numberOfElements1;j++){
 297         elementList.insert(number1[j]) ;
 298       }
 299
 300       for(int j=0;j<numberOfElements2;j++){
 301         elementList.insert(number2[j]) ;
 302       }
 303
 304       //create the array !
 305       int newNumberOfElements = elementList.size() ;
 306
 307       SCRUTE(newNumberOfElements);
 308
 309       numberOfElements[i] = newNumberOfElements ;
 310       int * tmp_arrayNew = new int[newNumberOfElements];
 311
 312       set<int>::iterator its ;
 313       for(its=elementList.begin();its!=elementList.end(); its++) {
 314         tmp_arrayNew[ii]=*its ;
 315         ii++;
 316       }
 317
 318       delete[] tmp_array[i] ;
 319       tmp_array[i] = tmp_arrayNew ;
 320       _totalNumberOfElements-=(numberOfElementsInType-newNumberOfElements) ;
 321
 322     } else
 323       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"ERROR")) ;
 324   }
 325   delete[] whereIsType ;
 326   delete[] tmp_NumberOfElementsInType ;
 327   delete [] _numberOfElements;
 328
 329   _numberOfElements = numberOfElements;
 330
 331   _numberOfGeometricType = it ;
 332   medGeometryElement * geometricType=_geometricType ;
 333   _geometricType = new medGeometryElement[it] ;
 334   int * numberOfGaussPoint=_numberOfGaussPoint ;
 335   _numberOfGaussPoint= new int[it] ;
 336 //    int * geometricTypeNumber=_geometricTypeNumber ;
 337 //    _geometricTypeNumber = new int[it] ;
 338
 339 //    MEDSKYLINEARRAY* numberNew = new MEDSKYLINEARRAY(it,_totalNumberOfElements);
 340 //    int * numberIndex = numberNew->getIndex() ;
 341   int size = _mesh->getNumberOfElements(_entity,MED_ALL_ELEMENTS);
 342   if (_totalNumberOfElements == size) _isOnAllElts = true;
 343
 344   int * numberValue = new int[_totalNumberOfElements] ;
 345   int * numberIndex = new int[it+1] ;
 346   numberIndex[0]=1;
 347   for (int i=0;i<it;i++) {
 348     memcpy(numberValue+numberIndex[i]-1,tmp_array[i],sizeof(int)*_numberOfElements[i]) ;
 349     delete[] tmp_array[i] ;
 350     numberIndex[i+1]=numberIndex[i]+_numberOfElements[i] ;
 351
 352     _geometricType[i]=myType[i] ;
 353     _numberOfGaussPoint[i]=1 ;
 354   }
 355   if ( _number != (MEDSKYLINEARRAY *) NULL) delete _number ;
 356   //_number = numberNew ;
 357   _number = new MEDSKYLINEARRAY(it,_totalNumberOfElements,numberIndex,numberValue);
 358   delete[] numberIndex;
 359
 360   delete[] numberValue;
 361
 362   delete[] myType ;
 363   delete[] tmp_array ;
 364
 365   delete[] geometricType ;
 366   delete[] numberOfGaussPoint ;
 367 //    delete[] geometricTypeNumber ;
 368 //  delete[] numberOfElements ;
 369
 370
 371
 372
 373   MESSAGE(LOC<<"Printing of the object SUPPORT blended "<< *this);
 374
 375
 376
 377
 378   END_OF(LOC);
 379 };
 380
 381 /*!
 382     This function allows the user to set a support not on all entities Entity,
 383     it should be used after an initialisation with the constructor
 384     SUPPORT(MESH* Mesh, string Name="", medEntityMesh Entity=MED_CELL) and
 385     after the call to the function setAll(false).
 386     It allocates and initialises all the attributs of the class SUPPORT.
 387  */
 388
 389 //-------------------
 390 void SUPPORT::setpartial(string Description, int NumberOfGeometricType,
 391                          int TotalNumberOfElements,
 392                          medGeometryElement *GeometricType,
 393                          int *NumberOfElements, int *NumberValue)
 394 //-------------------
 395 {
 396   const char * LOC = "SUPPORT::setpartial(string , int , int , medGeometryElement * , int * , int *) : " ;
 397   BEGIN_OF(LOC) ;
 398
 399   _isOnAllElts = false ;
 400
 401   _description=Description;
 402
 403   _numberOfGeometricType=NumberOfGeometricType;
 404
 405   if (_geometricType!=NULL) delete[] _geometricType ;
 406   _geometricType = new medGeometryElement[NumberOfGeometricType];
 407   if (_numberOfElements!=NULL) delete[] _numberOfElements ;
 408   _numberOfElements = new int[NumberOfGeometricType];
 409   _totalNumberOfElements = TotalNumberOfElements;
 410   if (_numberOfGaussPoint!=NULL) delete[] _numberOfGaussPoint ;
 411   _numberOfGaussPoint = new int[NumberOfGeometricType];
 412   int * index = new int[_numberOfGeometricType+1];
 413   index[0]=1;
 414   for (int i=0;i<_numberOfGeometricType;i++) {
 415     _geometricType[i] = GeometricType[i] ;
 416     _numberOfElements[i] = NumberOfElements[i] ;
 417     _numberOfGaussPoint[i] = 1 ;
 418     index[i+1] = index[i]+NumberOfElements[i] ;
 419   }
 420
 421   if (_number!=NULL) delete _number ;
 422   _number = new MEDSKYLINEARRAY(_numberOfGeometricType,_totalNumberOfElements,index,NumberValue);
 423
 424   delete[] index ;
 425
 426   END_OF(LOC);
 427 };
 428
 429
 430 /*!
 431   This method gets the boundary elements of the mesh. The support has to be
 432   build using the constructor SUPPORT(MESH *,string, medEntityMesh) or
 433   SUPPORT() followed by setMesh(MESH*) setName(string) and
 434   setEntity(medEntityMesh) before using this method.
 435 */
 436 //-------------------
 437 void SUPPORT::getBoundaryElements() throw (MEDEXCEPTION)
 438 //-------------------
 439 {
 440   const char * LOC = "SUPPORT::getBoundaryElements() : " ;
 441   BEGIN_OF(LOC) ;
 442
 443   if (_mesh == (MESH*)NULL) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"You shlould use the appropriate SUPPORT Constructor before calling this method"));
 444
 445   int spaceDimension = _mesh->getSpaceDimension();
 446
 447   if (spaceDimension == 3)
 448     if (_entity != MED_FACE)
 449       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<_entity<<" !"));
 450   if (spaceDimension == 2)
 451     if (_entity != MED_EDGE)
 452       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<_entity<<" !"));
 453
 454   setAll(false);
 455
 456   const int * myConnectivityValue = _mesh->getReverseConnectivity(MED_DESCENDING) ;
 457   const int * myConnectivityIndex = _mesh->getReverseConnectivityIndex(MED_DESCENDING) ;
 458   int numberOf = _mesh->getNumberOfElements(_entity,MED_ALL_ELEMENTS) ;
 459   list<int> myElementsList ;
 460   int size = 0 ;
 461   SCRUTE(numberOf) ;
 462   for (int i=0 ; i<numberOf; i++)
 463     if (myConnectivityValue[myConnectivityIndex[i]] == 0) {
 464       SCRUTE(i+1) ;
 465       myElementsList.push_back(i+1) ;
 466       size++ ;
 467     }
 468   SCRUTE(size) ;
 469   // Well, we must know how many geometric type we have found
 470   int * myListArray = new int[size] ;
 471   int id = 0 ;
 472   list<int>::iterator myElementsListIt ;
 473   for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
 474     myListArray[id]=(*myElementsListIt) ;
 475     SCRUTE(id);
 476     SCRUTE(myListArray[id]);
 477     id ++ ;
 478   }
 479
 480   int numberOfGeometricType ;
 481   medGeometryElement* geometricType ;
 482   int * numberOfGaussPoint ;
 483   int * geometricTypeNumber ;
 484   int * numberOfElements ;
 485   //MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY() ;
 486   int * mySkyLineArrayIndex ;
 487
 488   int numberOfType = _mesh->getNumberOfTypes(_entity) ;
 489   if (numberOfType == 1) { // wonderfull : it's easy !
 490     numberOfGeometricType = 1 ;
 491     geometricType = new medGeometryElement[1] ;
 492     const medGeometryElement *  allType = _mesh->getTypes(_entity);
 493     geometricType[0] = allType[0] ;
 494     numberOfGaussPoint = new int[1] ;
 495     numberOfGaussPoint[0] = 1 ;
 496     geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
 497     geometricTypeNumber[0] = 0 ;
 498     numberOfElements = new int[1] ;
 499     numberOfElements[0] = size ;
 500     mySkyLineArrayIndex = new int[2] ;
 501     mySkyLineArrayIndex[0]=1 ;
 502     mySkyLineArrayIndex[1]=1+size ;
 503   }
 504   else {// hemmm
 505     map<medGeometryElement,int> theType ;
 506     for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
 507       medGeometryElement myType = _mesh->getElementType(_entity,*myElementsListIt) ;
 508       if (theType.find(myType) != theType.end() )
 509         theType[myType]+=1 ;
 510       else
 511         theType[myType]=1 ;
 512     }
 513     numberOfGeometricType = theType.size() ;
 514     geometricType = new medGeometryElement[numberOfGeometricType] ;
 515     //const medGeometryElement *  allType = _mesh->getTypes(_entity); !! UNUSED VARIABLE !!
 516     numberOfGaussPoint = new int[numberOfGeometricType] ;
 517     geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
 518     numberOfElements = new int[numberOfGeometricType] ;
 519     mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
 520     int index = 0 ;
 521     mySkyLineArrayIndex[0]=1 ;
 522     map<medGeometryElement,int>::iterator theTypeIt ;
 523     for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
 524       geometricType[index] = (*theTypeIt).first ;
 525       numberOfGaussPoint[index] = 1 ;
 526       geometricTypeNumber[index] = 0 ;
 527       numberOfElements[index] = (*theTypeIt).second ;
 528       mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfElements[index] ;
 529       index++ ;
 530     }
 531   }
 532   //mySkyLineArray->setMEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
 533   MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
 534
 535   setNumberOfGeometricType(numberOfGeometricType) ;
 536   //  setGeometricType(geometricType) ;
 537   //  setNumberOfGaussPoint(numberOfGaussPoint) ;
 538   for (int i=0;i<numberOfGeometricType;i++)
 539     {
 540       _numberOfGaussPoint[i] = numberOfGaussPoint[i];
 541       _geometricType[i] = geometricType[i];
 542     }
 543
 544   setNumberOfElements(numberOfElements) ;
 545   setTotalNumberOfElements(size) ;
 546   //  setNumber(mySkyLineArray) ;
 547
 548   _number = new MEDSKYLINEARRAY(numberOfGeometricType,size);
 549
 550   _number->setIndex(mySkyLineArrayIndex);
 551
 552   for (int i=0;i<size;i++)
 553     {
 554       _number->setIndexValue(i+1,myListArray[i]);
 555     }
 556
 557   delete[] numberOfElements;
 558   delete[] geometricTypeNumber;
 559   delete[] numberOfGaussPoint;
 560   delete[] geometricType;
 561   delete[] mySkyLineArrayIndex;
 562   delete[] myListArray;
 563   delete mySkyLineArray;
 564
 565   END_OF(LOC) ;
 566 }
 567
 568 /*!
 569   intersect the given SUPPORT mySupport into the calling SUPPORT object.
 570 */
 571 //-------------------
 572 void SUPPORT::intersecting(SUPPORT * mySupport) throw (MEDEXCEPTION)
 573 //-------------------
 574 {
 575   const char * LOC = "SUPPORT::intersecting(SUPPORT *) : " ;
 576   BEGIN_OF(LOC) ;
 577
 578   MESSAGE(LOC<< "SUPPORT entry : " << *mySupport) ;
 579
 580   MESSAGE(LOC<< "SUPPORT (calling object) : " << *this) ;
 581
 582   // on same entity :
 583   if ( _entity != mySupport->getEntity() )
 584     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Entities are different !"));
 585
 586   int * tmp_NumberOfElementsInType = new int[MED_NBR_GEOMETRIE_MAILLE];
 587   medGeometryElement * myType = new medGeometryElement[MED_NBR_GEOMETRIE_MAILLE];
 588   int * whereIsType = new int[MED_NBR_GEOMETRIE_MAILLE];
 589   //MESH_ENTITIES myMeshEntities() ;
 590   list<MED_FR::med_geometrie_element>::const_iterator listIt ;
 591   int it=0 ;
 592   for(listIt=(MED_FR::meshEntities[(MED_FR::med_entite_maillage)_entity]).begin();listIt!=(MED_FR::meshEntities[(MED_FR::med_entite_maillage)_entity]).end();listIt++) {
 593     tmp_NumberOfElementsInType[it]=0;
 594     whereIsType[it]=0 ;
 595     myType[it]= MED_NONE;
 596     try {
 597       tmp_NumberOfElementsInType[it]+=getNumberOfElements((medGeometryElement)(*listIt)) ;
 598       whereIsType[it]+=1 ;
 599     }
 600     catch (const MEDEXCEPTION & ex) {};
 601     try {
 602       tmp_NumberOfElementsInType[it]+=mySupport->getNumberOfElements((medGeometryElement)(*listIt)) ;
 603       whereIsType[it]+=2 ;
 604     }
 605     catch (const MEDEXCEPTION & ex) {};
 606     if (whereIsType[it]==3) {
 607       myType[it]=(medGeometryElement)(*listIt) ;
 608       it++;
 609     }
 610   }
 611
 612   MESSAGE("it = "<< it);
 613
 614   // set new value :
 615   int * numberOfElements=_numberOfElements ;
 616   _numberOfElements = new int[it] ;
 617   _totalNumberOfElements = 0 ;
 618   //int totalSize = 0 ;
 619   int ** tmp_array = new (int*)[it];
 620   for (int i=0;i<it;i++) {
 621     int numberOfElementsInType = tmp_NumberOfElementsInType[i] ;
 622     _numberOfElements[i] = numberOfElementsInType ;
 623     tmp_array[i] = new int[numberOfElementsInType] ;
 624     _totalNumberOfElements+=numberOfElementsInType ;
 625     if (whereIsType[i] == 3) {
 626       const int * number1 = getNumber(myType[i]) ;
 627       const int * number2 = mySupport->getNumber(myType[i]) ;
 628
 629       SCRUTE(number1);
 630       SCRUTE(number2);
 631
 632       int numberOfElements1 = numberOfElements[i] ;
 633       int numberOfElements2 = mySupport->getNumberOfElements(myType[i]) ;
 634
 635       SCRUTE(numberOfElements1);
 636       SCRUTE(numberOfElements2);
 637
 638       set<int> setList1(number1,number1+numberOfElements1);
 639       set<int> setList2(number2,number2+numberOfElements2);
 640
 641       for(set<int>::iterator its=setList1.begin();its!=setList1.end(); its++)
 642         {
 643           MESSAGE("Number1 " << *its);
 644         }
 645
 646       for(set<int>::iterator its=setList2.begin();its!=setList2.end(); its++)
 647         {
 648           MESSAGE("Number2 " << *its);
 649         }
 650
 651       set<int> setListIntersect;
 652
 653       set_intersection(setList1.begin(),setList1.end(),setList2.begin(),
 654                        setList2.end(),inserter(setListIntersect,
 655                                                setListIntersect.begin()));
 656
 657       for(set<int>::iterator its=setListIntersect.begin();
 658           its!=setListIntersect.end(); its++)
 659         {
 660           MESSAGE("Number1 intersect Number2 " << *its);
 661         }
 662
 663       int newNumberOfElements = setListIntersect.size() ;
 664
 665       SCRUTE(newNumberOfElements);
 666
 667       _numberOfElements[i] = newNumberOfElements ;
 668       int * tmp_arrayNew = new int[newNumberOfElements];
 669
 670       int ii = 0 ;
 671
 672       for(set<int>::iterator its=setListIntersect.begin();
 673           its!=setListIntersect.end(); its++) {
 674         tmp_arrayNew[ii]=*its ;
 675         SCRUTE(tmp_arrayNew[ii]);
 676         ii++;
 677       }
 678
 679       delete[] tmp_array[i] ;
 680       tmp_array[i] = tmp_arrayNew ;
 681       _totalNumberOfElements-=(numberOfElementsInType-newNumberOfElements) ;
 682
 683     } else
 684       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"ERROR")) ;
 685   }
 686   delete[] whereIsType ;
 687   delete[] tmp_NumberOfElementsInType ;
 688
 689   _numberOfGeometricType = it ;
 690   medGeometryElement * geometricType=_geometricType ;
 691   _geometricType = new medGeometryElement[it] ;
 692   int * numberOfGaussPoint=_numberOfGaussPoint ;
 693   _numberOfGaussPoint= new int[it] ;
 694
 695   int size = _mesh->getNumberOfElements(_entity,MED_ALL_ELEMENTS);
 696   if (_totalNumberOfElements == size) _isOnAllElts = true;
 697
 698   int * numberValue = new int[_totalNumberOfElements] ;
 699   int * numberIndex = new int[it+1] ;
 700   numberIndex[0]=1;
 701   for (int i=0;i<it;i++) {
 702     memcpy(numberValue+numberIndex[i]-1,tmp_array[i],sizeof(int)*_numberOfElements[i]) ;
 703     delete[] tmp_array[i] ;
 704     numberIndex[i+1]=numberIndex[i]+_numberOfElements[i] ;
 705
 706     _geometricType[i]=myType[i] ;
 707     _numberOfGaussPoint[i]=1 ;
 708   }
 709   if ( _number != (MEDSKYLINEARRAY *) NULL) delete _number ;
 710
 711   _number = new MEDSKYLINEARRAY(it,_totalNumberOfElements,numberIndex,numberValue);
 712   delete[] numberIndex;
 713
 714   delete[] numberValue;
 715
 716   delete[] myType ;
 717   delete[] tmp_array ;
 718
 719   delete[] geometricType ;
 720   delete[] numberOfGaussPoint ;
 721 //    delete[] geometricTypeNumber ;
 722   delete[] numberOfElements ;
 723
 724   END_OF(LOC);
 725 };