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