src/MEDMEM/MEDMEM_Support.cxx

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