src/MEDMEM/MEDMEM_Support.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 /*
  24  File Support.cxx
  25 */
  26
  27 #include "MEDMEM_Support.hxx"
  28 #include "MEDMEM_DriversDef.hxx"
  29 #include "MEDMEM_GMesh.hxx"
  30 #include "MEDMEM_Meshing.hxx"
  31
  32 #include <set>
  33 #include <algorithm>
  34 #include <list>
  35
  36 using namespace std;
  37 using namespace MED_EN;
  38 using namespace MEDMEM;
  39
  40 #define MED_NBR_GEOMETRIE_MAILLE 15
  41
  42 /*!
  43 \defgroup SUPPORT_general General information
  44
  45 \defgroup SUPPORT_creation Creation methods
  46 The creation of a support requires a number of information
  47 which is supplied to the MedMem library with the following methods.
  48 When the support is defined on all elements, the creation method is
  49 very simple, for the element list is implicitly defined.
  50
  51 \defgroup SUPPORT_query Query methods
  52
  53 \defgroup SUPPORT_constructors Constructors
  54
  55 \defgroup SUPPORT_advanced Advanced methods
  56
  57 */
  58
  59
  60
  61 /* This class is a generic class for family and group */
  62
  63 /*!
  64   Constructor.
  65 */
  66 //--------------------------------------------------------------------------
  67 SUPPORT::SUPPORT(): _name(""),  _description("None"), _mesh((GMESH*)NULL),
  68                     _entity(MED_CELL), _numberOfGeometricType(0),
  69                     _isOnAllElts(false),
  70                     _totalNumberOfElements(0),
  71                     _number((MEDSKYLINEARRAY*)NULL),
  72                     _number_fromfile(0)
  73   //--------------------------------------------------------------------------
  74 {
  75   MESSAGE_MED("SUPPORT::SUPPORT()");
  76 }
  77
  78 /*!
  79 \addtogroup SUPPORT_constructors
  80 @{
  81 */
  82
  83 /*!
  84   Constructor of a support lying on mesh \a Mesh. By default,
  85   the support lies on all elements of type \a Entity.
  86   Partial support can be described using \a setpartial method.
  87
  88 \param Mesh Pointer to the mesh on which the support lies
  89 \param Name Support name (should not exceed MED_TAILLE_NOM as defined in Med - i.e. 32 characters)
  90 \param Entity Entity type of the support (MED_CELL,MED_FACE,MED_EDGE, MED_NODE)
  91 */
  92 //--------------------------------------------------------------------------
  93 // SUPPORT::SUPPORT(MESH* Mesh, string Name/*=""*/, medEntityMesh Entity/*=MED_CELL*/):
  94 //   _name(Name), _description("None"), _mesh(Mesh), _entity(Entity),
  95 //   _numberOfGeometricType(0), _isOnAllElts(true),
  96 //   _totalNumberOfElements(0), _number((MEDSKYLINEARRAY*)NULL),_number_fromfile(0)
  97 //   //--------------------------------------------------------------------------
  98 // {
  99 //   MESSAGE_MED("SUPPORT::SUPPORT(MESH*Mesh,string Name,medEntityMesh Entity)");
 100 //   if(_mesh)
 101 //     _mesh->addReference();
 102 //   update() ;
 103 // }
 104
 105 /*!
 106   Copy constructor.
 107 */
 108 //--------------------------------------------------------------------------
 109 SUPPORT::SUPPORT(const SUPPORT & m):_number_fromfile(0)
 110   //--------------------------------------------------------------------------
 111 {
 112   const char* LOC = "SUPPORT::SUPPORT(SUPPORT & m) : ";
 113   BEGIN_OF_MED(LOC);
 114
 115   _name = m._name ;
 116   _description = m._description ;
 117   _mesh = m._mesh ; // on recopie uniquement l'adresse
 118   if(_mesh)
 119     _mesh->addReference();
 120   _meshName = m._meshName;
 121   _entity = m._entity;
 122   _numberOfGeometricType = m._numberOfGeometricType;
 123
 124   if (m._geometricType)
 125     _geometricType.set(_numberOfGeometricType,m._geometricType);
 126
 127   _isOnAllElts = m._isOnAllElts;
 128
 129   if (m._numberOfElements)
 130     _numberOfElements.set(_numberOfGeometricType,m._numberOfElements);
 131
 132   _totalNumberOfElements = m._totalNumberOfElements;
 133
 134   if ( m._number ) // m may be not filled SUPPORTClient
 135     _number = new MEDSKYLINEARRAY(* m._number);
 136   else
 137     _number = (MEDSKYLINEARRAY *) NULL;
 138
 139   _profilNames=m._profilNames;
 140
 141   END_OF_MED(LOC);
 142 }
 143
 144 /*!
 145   @}
 146 */
 147
 148 /*!\ifnot MEDMEM_ug
 149   Affectation operator. operator = perform et deep copy except for attribute _mesh
 150 \endif
 151 */
 152 //--------------------------------------------------------------------------
 153 SUPPORT & SUPPORT::operator=(const SUPPORT & m)
 154   //--------------------------------------------------------------------------
 155 {
 156   const char* LOC = "SUPPORT::operator=(const SUPPORT & m) : ";
 157   BEGIN_OF_MED(LOC);
 158
 159   if ( this == &m ) return *this;
 160
 161   _name = m._name;
 162   _description = m._description;
 163   if(m._mesh!=_mesh)//setMesh not used here due to _meshName update is this...
 164     {
 165       if(_mesh)
 166         _mesh->removeReference();
 167       _mesh=m._mesh;
 168       if(_mesh)
 169         _mesh->addReference();
 170     }
 171   _entity = m._entity;
 172   _numberOfGeometricType = m._numberOfGeometricType;
 173   if (m._geometricType)
 174     _geometricType.set(_numberOfGeometricType,m._geometricType);
 175   _isOnAllElts = m._isOnAllElts;
 176   if (m._numberOfElements)
 177     _numberOfElements.set(_numberOfGeometricType,m._numberOfElements);
 178   _totalNumberOfElements = m._totalNumberOfElements;
 179
 180   if (_number) delete _number;
 181   if  ( m._number ) // m may be not filled SUPPORTClient
 182     _number = new MEDSKYLINEARRAY(* m._number);
 183   else
 184     _number = (MEDSKYLINEARRAY *) NULL;
 185
 186   _profilNames=m._profilNames;
 187
 188   END_OF_MED(LOC);
 189   return *this;
 190 }
 191
 192 /*!
 193   Destructor.
 194 */
 195 //-----------------
 196 SUPPORT::~SUPPORT()
 197   //-----------------
 198 {
 199   MESSAGE_MED("Destructeur ~SUPPORT()");
 200   clearDataOnNumbers();
 201   if(_mesh)
 202     _mesh->removeReference();
 203 }
 204
 205 /*!
 206   operator <<.
 207 */
 208 //--------------------------------------------------
 209 ostream & MEDMEM::operator<<(ostream &os, const SUPPORT &my)
 210   //--------------------------------------------------
 211 {
 212   os << "Name : "<< my.getName() << endl ;
 213   os << "Description : "<< my.getDescription() << endl ;
 214   os << "Mesh ptr : ";
 215   if (my.getMesh() == NULL)
 216     os << " Mesh not defined." << endl ;
 217   else
 218     os << " Mesh defined." << endl;
 219   os << "MeshName : ";
 220   os << my.getMeshName() << endl ;
 221   os << "Entity : "<<entNames[my._entity] << endl;
 222   os << "Entity list : "<< endl;
 223   if ( my._isOnAllElts )
 224     os << "Is on all entities."<< endl;
 225   else {
 226     os << "Is not on all entities. "<< endl;
 227     if ( my._number )  // m may be not filled SUPPORTClient
 228       os << *my.getNumber(MED_ALL_ELEMENTS);
 229   }
 230   int numberoftypes = my._numberOfGeometricType ;
 231   os << "NumberOfTypes : "<<numberoftypes<<endl;
 232   PointerOf<medGeometryElement> types = my._geometricType;
 233   for (int j=0;j<numberoftypes;j++) {
 234     int numberOfElements = my._numberOfElements ? my._numberOfElements[j] : -1;
 235     os << "    On Type "<<geoNames[types[j]]
 236        <<" : there is(are) "<<numberOfElements<<" element(s) and " <<endl;
 237   }
 238   int nbProfilNames = my._profilNames.size();
 239   os<<"Number of profil names = "<<nbProfilNames<<endl;
 240   for(int j=0; j<nbProfilNames; j++) {
 241     os<<"    Profil Name N"<<j+1<<" = "<<my._profilNames[j]<<endl;
 242   }
 243   return os ;
 244 }
 245
 246 /*!
 247   Updade the SUPPORT attributs with right MESH information.
 248
 249   It has an effect only if SUPPORT is on all elements.
 250
 251   No more need in future release.
 252 */
 253 //-------------------
 254 void SUPPORT::update()
 255 //-------------------
 256 {
 257   const char* LOC = "SUPPORT::update() : ";
 258   BEGIN_OF_MED(LOC);
 259
 260   if (_isOnAllElts && _mesh)
 261     {
 262       if (_entity == MED_NODE)
 263         {
 264           // BEGIN Issue 0020804: [CEA 399] Memory corruption ... in MEDMEMCppTest
 265           //_numberOfGeometricType = 1;
 266           setNumberOfGeometricType(1);
 267           // END Issue 0020804
 268
 269           // BEGIN Issue 0020633: [CEA] Pb with 3D field creation fron another
 270           // Use setGeometricType() in order to get _profilNames updated
 271           //_geometricType.set(1);
 272           //_geometricType[0]=MED_POINT1;
 273           const MED_EN::medGeometryElement type = MED_NONE;
 274           setGeometricType( & type );
 275           // END Issue 0020633: [CEA] Pb with 3D field creation fron another
 276           _numberOfElements.set(1);
 277           _numberOfElements[0]=_mesh->getNumberOfNodes();
 278           _totalNumberOfElements=_numberOfElements[0];
 279         }
 280       else
 281         { // we duplicate information from _mesh
 282           // BEGIN Issue 0020804: [CEA 399] Memory corruption ... in MEDMEMCppTest
 283           setNumberOfGeometricType(_mesh->getNumberOfTypes(_entity));
 284           // END Issue 0020804
 285           // BEGIN Issue 0020633: [CEA] Pb with 3D field creation fron another
 286           if ( const medGeometryElement *  allType = _mesh->getTypes(_entity))
 287             setGeometricType( allType );
 288           // END Issue 0020633: [CEA] Pb with 3D field creation fron another
 289           _numberOfElements.set(_numberOfGeometricType);
 290           _totalNumberOfElements=0;
 291           for (int i=0;i<_numberOfGeometricType;i++)
 292             {
 293               _numberOfElements[i]=_mesh->getNumberOfElements(_entity,_geometricType[i]) ;
 294               _totalNumberOfElements+=_numberOfElements[i];
 295             }
 296         }
 297
 298       if (_totalNumberOfElements <= 0)
 299         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"We have found no element for this support !"));
 300       // set _number (issue 0021167)
 301       {
 302         vector<int> nums( _totalNumberOfElements );
 303         for ( unsigned i = 0; i < nums.size(); ++i )
 304           nums[i] = i+1;
 305
 306         vector<int> index( _numberOfGeometricType + 1 );
 307         index[0] = 1;
 308         for ( int i = 0; i < _numberOfGeometricType; ++i )
 309           index[i+1] = index[i] + _numberOfElements[i];
 310
 311         setNumber( & index[0], & nums[0] );
 312       }
 313     }
 314   END_OF_MED(LOC);
 315 }
 316
 317 /*!
 318   Get the field value index (in fortran mode) from the support global number.
 319   Becareful, it doesn't take care of the field number of components
 320 */
 321 //-------------------
 322 int SUPPORT::getValIndFromGlobalNumber(const int number) const throw (MEDEXCEPTION)
 323 //-------------------
 324 {
 325   const char * LOC="getValIndFromGlobalNumber(const int number) : ";
 326
 327   if (_isOnAllElts) return number;
 328
 329   int nbOfEltsThis    = getNumberOfElements(MED_ALL_ELEMENTS);
 330
 331   const int *eltsThis = _number->getValue();
 332
 333   int iThis;
 334   bool found=false;
 335
 336   for(iThis=0;iThis<nbOfEltsThis && !found;)
 337     if(eltsThis[iThis]==number)
 338     {
 339       found = true;
 340       int valInd = iThis+1;
 341       return valInd;
 342     }
 343     else
 344       iThis++;
 345
 346   if(!found)
 347     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't find the global number |"
 348                                  << number << "| in Support |"
 349                                  << getName() << "|" ));
 350
 351   // It should never arrive here !!
 352   return 0;
 353 }
 354
 355 /*!
 356 \addtogroup SUPPORT_advanced
 357 @{
 358
 359 */
 360
 361 /*!
 362   Blends the given SUPPORT mySupport into the calling object SUPPORT.
 363   Example :
 364 \verbatim
 365 SUPPORT mySupport ;
 366 SUPPORT myOtherSupport ;
 367 ...
 368 mySupport.blending(myOtherSupport) ;
 369 \endverbatim
 370 Support \a mySupport now contains a union of the elements originally
 371 contained in \a mySupport and \a myOtherSupport.
 372 */
 373 //-------------------
 374 void SUPPORT::blending(const SUPPORT * mySupport) throw (MEDEXCEPTION)
 375   //-------------------
 376 {
 377   const char * LOC="SUPPORT::blending(SUPPORT *) : ";
 378   BEGIN_OF_MED(LOC);
 379   if (_entity!=mySupport->getEntity())
 380     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Entities are different !"));
 381   if(!(*_mesh == *mySupport->getMesh()))
 382     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Mesh are different !"));
 383   if(_isOnAllElts)
 384     return;
 385   if(mySupport->isOnAllElements())
 386   {
 387     *this=*mySupport;
 388     return;
 389   }
 390   if(mySupport->_totalNumberOfElements==0)
 391     return;
 392   const int *ids=getNumber(MED_ALL_ELEMENTS);
 393   set<int> idsSet(ids,ids+getNumberOfElements(MED_ALL_ELEMENTS));
 394   const int *idsMySupport=mySupport->getNumber(MED_ALL_ELEMENTS);
 395   int mySupportSize=mySupport->getNumberOfElements(MED_ALL_ELEMENTS);
 396   set<int>::iterator iter;
 397   for(int i=0;i<mySupportSize;i++)
 398     idsSet.insert(idsMySupport[i]);
 399   int size=idsSet.size();
 400
 401   if(size!=0)
 402   {
 403     list<int> idsList;
 404     for(iter=idsSet.begin();iter!=idsSet.end();iter++)
 405       idsList.push_back(*iter);
 406
 407     MESSAGE_MED(LOC << " size Set " << idsSet.size() << " size List " << idsList.size());
 408
 409     if(_entity==MED_NODE)
 410       fillFromNodeList(idsList);
 411     else
 412       fillFromElementList(idsList);
 413   }
 414   else
 415     clearDataOnNumbers();
 416   END_OF_MED(LOC);
 417 }
 418 /*!  @}  */
 419
 420 /*!
 421 \addtogroup SUPPORT_creation
 422 @{
 423 */
 424
 425 /*!
 426   This function allows the user to set a support not on all entities Entity,
 427   it should be used after setting the mesh and the entity attributes at least.
 428
 429 \param Description string describing the support for information purposes (should not exceed MED_TAILLE_DESC length - i.e. 200 characters)
 430 \param NumberOfGeometricType number of geometric types contained in the support
 431 \param TotalNumberOfElements number of elements in the support
 432 \param GeometricType array describing the geometric types (must be consistent with the entity that was passed as an argument to the support constructor)
 433 \param NumberOfElements array describing the number of elements for each type
 434 \param NumberValue array of IDs of the elements that constitute the group.
 435
 436 The following example refers to the mesh given in the mesh connectivity example.
 437 It creates a group containing the two cells on the right (the quadratic triangle and the quadrangle on the right).
 438
 439 \verbatim
 440 // creating SUPPORT on cells with one value per cell
 441 right_group = new SUPPORT;
 442 right_group->setMesh(mesh);
 443 right_group->setEntity( MED_CELL );
 444 right_group->setName("right group");
 445
 446 string description = "partial support";
 447 int number_of_types=2;
 448 int number_of_elements=2;
 449 medGeometryElement geom_types[2]={MED_QUAD4, MED_TRIA6};
 450 int number_of_elem_per_type[2]={1,1};
 451 int number_value[2]={3,4};
 452
 453 //defining the region of the support
 454 right_group->setpartial(description, number_of_types,
 455 number_of_elements, geom_types,
 456 number_of_elem_per_type, number_value);
 457 \endverbatim
 458
 459 When MED_POLYGON or MED_POLYHEDRA elements are included in the support,
 460 their global number should be given. For instance, on a mesh having ten MED_TRIA3
 461 and five MED_POLYGON, the number of the first polygonal element is 11.
 462 */
 463
 464 //-------------------
 465 void SUPPORT::setpartial(const std::string&         Description,
 466                          int                        NumberOfGeometricType,
 467                          int                        TotalNumberOfElements,
 468                          const medGeometryElement * GeometricType,
 469                          const int *                NumberOfElements,
 470                          const int *                NumberValue)
 471 //-------------------
 472 {
 473   const char * LOC = "SUPPORT::setpartial(string , int , int , medGeometryElement * , int * , int *) : " ;
 474   BEGIN_OF_MED(LOC) ;
 475
 476   _isOnAllElts = false ;
 477
 478   _description=Description;
 479
 480   _numberOfGeometricType=NumberOfGeometricType;
 481   _geometricType.set(NumberOfGeometricType);
 482   _numberOfElements.set(NumberOfGeometricType);
 483   _totalNumberOfElements = TotalNumberOfElements;
 484
 485   int * index = new int[_numberOfGeometricType+1];
 486   index[0]=1;
 487   int elemDim = -1;
 488   for (int i=0;i<_numberOfGeometricType;i++) {
 489     if(GeometricType[i]/100 != elemDim)
 490       {
 491         if(i==0)
 492           elemDim=GeometricType[i]/100;
 493         else if ( CELLMODEL_Map::retrieveCellModel( GeometricType[i] ).getDimension() !=
 494                   CELLMODEL_Map::retrieveCellModel( GeometricType[0] ).getDimension() )
 495           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"unhomogeneous geometric types (dimension) !"));
 496       }
 497     _geometricType[i] = GeometricType[i] ;
 498     _numberOfElements[i] = NumberOfElements[i] ;
 499     index[i+1] = index[i]+NumberOfElements[i] ;
 500   }
 501
 502   if (_number!=NULL) delete _number ;
 503   _number = new MEDSKYLINEARRAY(_numberOfGeometricType,_totalNumberOfElements,index,NumberValue);
 504
 505   delete[] index ;
 506
 507   // PAL16854(Partial support on nodes):
 508   // giving a default value to profile names
 509   vector<string> prof_names( NumberOfGeometricType);
 510   for (int itype=0; itype < NumberOfGeometricType; itype++)
 511   {
 512     ostringstream typestr;
 513     typestr<<_name<<"_type"<<_geometricType[itype];
 514     prof_names[itype]=typestr.str();
 515   }
 516   setProfilNames(prof_names);
 517
 518   END_OF_MED(LOC);
 519 }
 520
 521 /*! @}  */
 522
 523 /*!
 524 \ifnot MEDMEM_ug
 525 This function allows the user to set a support not on all entities Entity,
 526 it should be used after an initialisation of :
 527 SUPPORT(GMESH* Mesh, string Name="", medEntityMesh Entity=MED_CELL) and
 528 after calling  at least setGeometricType and perharps setEntity.
 529 It allocates and initialises all the attributs of the class SUPPORT but
 530 doesn't set a description, a SUPPORT name, a meshName and an associated GMESH.
 531 \endif
 532 */
 533
 534 //-------------------
 535 void SUPPORT::setpartial(MEDSKYLINEARRAY * number, bool shallowCopy) throw (MEDEXCEPTION)
 536   //-------------------
 537 {
 538   const char * LOC = "SUPPORT::setpartial(MEDSKYLINEARRAY * number) : " ;
 539   BEGIN_OF_MED(LOC) ;
 540
 541   if ( ! _geometricType )
 542     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"SUPPORT must contains"
 543                                  << " a geometric type list" )) ;
 544
 545   _numberOfGeometricType = number->getNumberOf();
 546
 547   _numberOfElements.set(_numberOfGeometricType);
 548
 549   for (int i=0; i< _numberOfGeometricType; i++)
 550     _numberOfElements[i] =  number->getNumberOfI(i+1);
 551
 552   _totalNumberOfElements = number->getLength();
 553
 554   _isOnAllElts = false ;
 555
 556   if (_number!=NULL) delete _number ;
 557
 558   if ( shallowCopy )
 559     _number = number;
 560   else
 561     _number = new MEDSKYLINEARRAY(*number);
 562
 563   END_OF_MED(LOC);
 564 }
 565
 566 void SUPPORT::setpartial_fromfile(MEDSKYLINEARRAY * number, bool shallowCopy) throw (MEDEXCEPTION)
 567   //-------------------
 568 {
 569   const char* LOC = "SUPPORT::setpartial_fromfile(MEDSKYLINEARRAY * number) : ";
 570   BEGIN_OF_MED(LOC);
 571
 572   if ( shallowCopy )
 573     _number_fromfile = number;
 574   else
 575     _number_fromfile = new MEDSKYLINEARRAY(*number);
 576
 577   END_OF_MED(LOC);
 578 }
 579
 580 void SUPPORT::setProfilNames(const std::vector<std::string>& profilNames) throw (MEDEXCEPTION){
 581
 582   const char * LOC = "SUPPORT::setProfilNames(vector<string> profilNames) : " ;
 583   BEGIN_OF_MED(LOC) ;
 584
 585   if ( _isOnAllElts )
 586     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"SUPPORT shouldn't be on all elements"
 587                                  << " while setting profil name list" )) ;
 588
 589   if ( ! _geometricType || _numberOfGeometricType==0 )
 590     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"SUPPORT must contains"
 591                                  << " a least one geometric type" )) ;
 592
 593   if ( ! _number )
 594     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"SUPPORT must contains"
 595                                  << " a profil number list before setting"
 596                                  << " the associated profil name list" )) ;
 597
 598   if ( ( (int)profilNames.size() != _number->getNumberOf() ) &&
 599        ( (int)profilNames.size() !=_numberOfGeometricType ) ) {
 600     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The profil name list size : "<< profilNames.size()
 601                                  << " must be equal to the number of geometric type : "
 602                                  <<  _numberOfGeometricType << " (_number->getNumberOf() : "
 603                                  << _number->getNumberOf() << " )"
 604                                  )) ;
 605
 606   }
 607
 608   _profilNames = profilNames;
 609
 610   END_OF_MED(LOC);
 611
 612 }
 613
 614 vector<string> SUPPORT::getProfilNames() const throw (MEDEXCEPTION)
 615 {
 616   return _profilNames;
 617 }
 618
 619 /*!
 620 \addtogroup SUPPORT_advanced
 621 @{
 622 */
 623
 624 /*!
 625   This method gets the boundary elements of the mesh. The support has to be
 626   build using SUPPORT() followed by setMesh(GMESH*) setName(string) and
 627   setEntity(medEntityMesh) before using this method.
 628 */
 629 //-------------------
 630 void SUPPORT::getBoundaryElements() throw (MEDEXCEPTION)
 631   //-------------------
 632 {
 633   const char * LOC = "SUPPORT::getBoundaryElements() : " ;
 634   BEGIN_OF_MED(LOC) ;
 635
 636   if (_mesh == (GMESH*)NULL) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"You shlould use the appropriate SUPPORT Constructor before calling this method"));
 637
 638   int spaceDimension = _mesh->getSpaceDimension();
 639
 640   if (spaceDimension == 3)
 641     if (_entity != MED_FACE)
 642       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<_entity<<" !"));
 643   if (spaceDimension == 2)
 644     if (_entity != MED_EDGE)
 645       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<_entity<<" !"));
 646
 647   setAll(false);
 648
 649   const MESH* mesh = _mesh->convertInMESH();
 650   const_cast<CONNECTIVITY*>
 651     (mesh->getConnectivityptr())->calculateFullDescendingConnectivity(MED_CELL);
 652   const int * myConnectivityValue = mesh->getReverseConnectivity(MED_DESCENDING) ;
 653   const int * myConnectivityIndex = mesh->getReverseConnectivityIndex(MED_DESCENDING) ;
 654   int numberOf = mesh->getNumberOfElements(_entity,MED_ALL_ELEMENTS) ;
 655   list<int> myElementsList ;
 656   int size = 0 ;
 657   SCRUTE_MED(numberOf) ;
 658   for (int i=0 ; i<numberOf; i++)
 659     if (myConnectivityValue[myConnectivityIndex[i]] == 0) {
 660       SCRUTE_MED(i+1) ;
 661       myElementsList.push_back(i+1) ;
 662       size++ ;
 663     }
 664   SCRUTE_MED(size) ;
 665   // Well, we must know how many geometric type we have found
 666   int * myListArray = new int[size] ;
 667   int id = 0 ;
 668   list<int>::iterator myElementsListIt ;
 669   for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
 670     myListArray[id]=(*myElementsListIt) ;
 671     SCRUTE_MED(id);
 672     SCRUTE_MED(myListArray[id]);
 673     id ++ ;
 674   }
 675
 676   int numberOfGeometricType ;
 677   medGeometryElement* geometricType ;
 678   int * geometricTypeNumber ;
 679   int * numberOfElements ;
 680   int * mySkyLineArrayIndex ;
 681
 682   int numberOfType = mesh->getNumberOfTypes(_entity) ;
 683   if (numberOfType == 1) { // wonderfull : it's easy !
 684     numberOfGeometricType = 1 ;
 685     geometricType = new medGeometryElement[1] ;
 686     const medGeometryElement *  allType = mesh->getTypes(_entity);
 687     geometricType[0] = allType[0] ;
 688     geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
 689     geometricTypeNumber[0] = 0 ;
 690     numberOfElements = new int[1] ;
 691     numberOfElements[0] = size ;
 692     mySkyLineArrayIndex = new int[2] ;
 693     mySkyLineArrayIndex[0]=1 ;
 694     mySkyLineArrayIndex[1]=1+size ;
 695   }
 696   else {// hemmm
 697     map<medGeometryElement,int> theType ;
 698     for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
 699       medGeometryElement myType = mesh->getElementType(_entity,*myElementsListIt) ;
 700       if (theType.find(myType) != theType.end() )
 701         theType[myType]+=1 ;
 702       else
 703         theType[myType]=1 ;
 704     }
 705     numberOfGeometricType = theType.size() ;
 706     geometricType = new medGeometryElement[numberOfGeometricType] ;
 707     geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
 708     numberOfElements = new int[numberOfGeometricType] ;
 709     mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
 710     int index = 0 ;
 711     mySkyLineArrayIndex[0]=1 ;
 712     map<medGeometryElement,int>::iterator theTypeIt ;
 713     for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
 714       geometricType[index] = (*theTypeIt).first ;
 715       geometricTypeNumber[index] = 0 ;
 716       numberOfElements[index] = (*theTypeIt).second ;
 717       mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfElements[index] ;
 718       index++ ;
 719     }
 720   }
 721   MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
 722
 723   setNumberOfGeometricType(numberOfGeometricType) ;
 724   setGeometricType(geometricType) ;
 725   setNumberOfElements(numberOfElements) ;
 726
 727   _number = new MEDSKYLINEARRAY(numberOfGeometricType,size);
 728
 729   _number->setIndex(mySkyLineArrayIndex);
 730
 731   for (int i=0;i<size;i++)
 732   {
 733     _number->setIndexValue(i+1,myListArray[i]);
 734   }
 735
 736   delete[] numberOfElements;
 737   delete[] geometricTypeNumber;
 738   delete[] geometricType;
 739   delete[] mySkyLineArrayIndex;
 740   delete[] myListArray;
 741   delete mySkyLineArray;
 742
 743   mesh->removeReference();
 744
 745   END_OF_MED(LOC);
 746 }
 747
 748 /*!
 749   Intersects \a mySupport into the calling SUPPORT object.
 750   If A.intersecting(B) is called, on output, \f$ A \f$ contains \f$A \cap B\f$.
 751 */
 752 //-------------------
 753 void SUPPORT::intersecting(const SUPPORT * mySupport) throw (MEDEXCEPTION)
 754 {
 755   const char * LOC="SUPPORT::intersecting(SUPPORT *) : ";
 756   BEGIN_OF_MED(LOC);
 757   if (_entity!=mySupport->getEntity())
 758     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Entities are different !"));
 759   if(!(*_mesh == *mySupport->getMesh()))
 760     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Mesh are different !"));
 761   if(mySupport->isOnAllElements())
 762     return;
 763   if(_isOnAllElts)
 764   {
 765     *this=*mySupport;
 766     return;
 767   }
 768   if(_totalNumberOfElements==0)
 769     return;
 770   const int *ids=getNumber(MED_ALL_ELEMENTS);
 771   set<int> idsSet(ids,ids+getNumberOfElements(MED_ALL_ELEMENTS));
 772   const int *idsMySupport=mySupport->getNumber(MED_ALL_ELEMENTS);
 773   int mySupportSize=mySupport->getNumberOfElements(MED_ALL_ELEMENTS);
 774   set<int> idsSetMySupport(idsMySupport,idsMySupport+mySupportSize);
 775   set<int>::iterator iter;
 776   list<int> idsList;
 777   for(iter=idsSet.begin();iter!=idsSet.end();iter++)
 778     if(idsSetMySupport.find(*iter)!=idsSetMySupport.end())
 779       idsList.push_back(*iter);
 780   int size=idsSet.size();
 781   int sizeList = idsList.size();
 782
 783   MESSAGE_MED(LOC << " size Set " << idsSet.size() << " size List " << idsList.size());
 784
 785   if(size!=0 && sizeList != 0)
 786   {
 787     if(_entity==MED_NODE)
 788       fillFromNodeList(idsList);
 789     else
 790       fillFromElementList(idsList);
 791   }
 792   else
 793   {
 794     clearDataOnNumbers();
 795   }
 796   END_OF_MED(LOC);
 797 }
 798 /*!  @}  */
 799
 800 /*!
 801   Method that cleans up all the fields related to _numbers. Defined for code factorization.
 802 */
 803 //--------------------------------------------------
 804 void MEDMEM::SUPPORT::clearDataOnNumbers()
 805   //--------------------------------------------------
 806 {
 807   _numberOfGeometricType=0;
 808   _totalNumberOfElements=0;
 809
 810   if(_number)
 811   {
 812     delete _number;
 813     _number=(MEDSKYLINEARRAY *) NULL;
 814   }
 815   if(_number_fromfile)
 816     {
 817       delete _number_fromfile;
 818       _number_fromfile=0;
 819     }
 820 }
 821
 822 /*!
 823   operator == This operator does not compare attributs _name and _description.
 824 */
 825 //--------------------------------------------------
 826 bool MEDMEM::SUPPORT::operator == (const SUPPORT &support) const
 827   //--------------------------------------------------
 828 {
 829
 830   const char* LOC = "bool SUPPORT::operator ==(const SUPPORT &support) const : ";
 831   BEGIN_OF_MED(LOC);
 832
 833   bool operatorReturn = false;
 834
 835   operatorReturn = (*_mesh == *support._mesh) && (_entity == support._entity) &&
 836     (_numberOfGeometricType == support._numberOfGeometricType) &&
 837     ((_isOnAllElts && support._isOnAllElts) || (!_isOnAllElts && !support._isOnAllElts)) &&
 838     (_totalNumberOfElements == support._totalNumberOfElements) &&
 839     (_profilNames.size() == support._profilNames.size());
 840
 841   if (operatorReturn)
 842   {
 843     if (!_isOnAllElts)
 844     {
 845       for (int i=0; i<_numberOfGeometricType; i++)
 846       {
 847         operatorReturn = operatorReturn &&
 848           (_geometricType[i] == support._geometricType[i]) &&
 849           (_numberOfElements[i] == support._numberOfElements[i]);
 850
 851         if (operatorReturn)
 852         {
 853           for (int j=0; j<_numberOfElements[i]; j++)
 854           {
 855             operatorReturn = operatorReturn &&
 856               (getNumber(_geometricType[i])[j] ==
 857                support.getNumber(_geometricType[i])[j]);
 858           }
 859         }
 860       }
 861     }
 862   }
 863
 864   END_OF_MED(LOC);
 865
 866   return operatorReturn;
 867 }
 868
 869 void SUPPORT::changeElementsNbs(medEntityMesh entity, const int *renumberingFromOldToNew)
 870 {
 871   if(entity != _entity)
 872     throw MEDEXCEPTION("SUPPORT::changeElementsNbs : Renumbering on a mismatch entity");
 873   list<int> newNbs;
 874   if(!_isOnAllElts)
 875   {
 876     const int *oldNbs=_number->getValue();
 877     for(int i=0;i<_totalNumberOfElements;i++)
 878       newNbs.push_back(renumberingFromOldToNew[oldNbs[i]-1]);
 879     newNbs.sort();
 880     fillFromElementList(newNbs);
 881   }
 882   else
 883     update();
 884 }
 885
 886 /*!
 887   operator == + in case false a test if coordinates and connectivity of _mesh and support->_mesh are the same
 888 */
 889 bool MEDMEM::SUPPORT::deepCompare(const SUPPORT &support) const
 890 {
 891   bool operatorReturn =(_entity == support._entity) &&
 892     (_numberOfGeometricType == support._numberOfGeometricType) &&
 893     ( (_isOnAllElts && support._isOnAllElts) || (!_isOnAllElts  && !support._isOnAllElts) ) &&
 894     (_totalNumberOfElements == support._totalNumberOfElements);
 895   if (operatorReturn)
 896   {
 897     if (!_isOnAllElts)
 898     {
 899       for (int i=0; i<_numberOfGeometricType && operatorReturn; i++)
 900       {
 901         operatorReturn = (_geometricType[i] == support._geometricType[i]) &&
 902           (_numberOfElements[i] == support._numberOfElements[i]);
 903         if (operatorReturn)
 904         {
 905           for (int j=0; j<_numberOfElements[i]; j++)
 906           {
 907             operatorReturn = (getNumber(_geometricType[i])[j] ==
 908                               support.getNumber(_geometricType[i])[j]);
 909           }
 910         }
 911       }
 912     }
 913   }
 914   if(operatorReturn)
 915     operatorReturn = ( bool(_mesh) == bool(support._mesh));
 916
 917   if(operatorReturn)
 918   {
 919     if(!(*_mesh == *support._mesh))
 920     {
 921       return _mesh->deepCompare(*support._mesh);
 922     }
 923   }
 924   return operatorReturn;
 925 }
 926
 927 /*!
 928   States if this is included in other.
 929 */
 930 bool MEDMEM::SUPPORT::belongsTo(const SUPPORT& other, bool deepCompare) const
 931 {
 932   if(!(*_mesh == *other._mesh))
 933   {
 934     if(!deepCompare)
 935       return false;
 936     if(!_mesh->deepCompare(*other._mesh))
 937       return false;
 938   }
 939   if(_entity!=other._entity)
 940     return false;
 941   if(other._isOnAllElts)
 942     return true;
 943   if(_isOnAllElts && !other._isOnAllElts)
 944     return false;
 945   if(_numberOfGeometricType>other._numberOfGeometricType)
 946     return false;
 947   for(int i=0; i<_numberOfGeometricType; i++)
 948   {
 949     medGeometryElement curGeomType=_geometricType[i];
 950     int iOther=-1;
 951     for(int j=0; j<other._numberOfGeometricType; j++)
 952       if(other._geometricType[j]==curGeomType)
 953         iOther=j;
 954     if(iOther==-1)
 955       return false;
 956     if(_numberOfElements[i]>other._numberOfElements[iOther])
 957       return false;
 958     const int *numbers1=_number->getI(i+1);
 959     const int *numbers2=other._number->getI(iOther+1);
 960     for (int k=0; k<_numberOfElements[i]; k++)
 961     {
 962       bool found=false;
 963       for(int l=0;l<other._numberOfElements[iOther] && !found;l++)
 964       {
 965         if(numbers1[k]==numbers2[l])
 966           found=true;
 967       }
 968       if(!found)
 969         return false;
 970     }
 971   }
 972   return true;
 973 }
 974 /*!
 975   Method used to sort array of id.
 976 */
 977 int compareId(const void *x, const void *y);
 978 int compareId(const void *x, const void *y)
 979 {
 980   const int *x1=(const int *)x;
 981   const int *y1=(const int *)y;
 982   if(*x1<*y1)
 983     return -1;
 984   else if(*x1>*y1)
 985     return 1;
 986   else
 987     return 0;
 988 }
 989
 990 /*!
 991   performs a common operation : Sub builds a sorted int array that is obtained by supression of all ids contained
 992   in array defined by (idsToSuppress,lgthIdsToSuppress) from array [start ... end]
 993   Example sub(0,7,{1,2,5},3) => {0,3,4,6,7} - WARNING returned list should be deallocated !
 994 */
 995 list<int> *MEDMEM::SUPPORT::sub(int start,int end,const int *idsToSuppress,int lgthIdsToSuppress)
 996 {
 997   int size=end-start+1;
 998   int sizeRet=size-lgthIdsToSuppress;
 999   list<int> *ret;
1000   if(sizeRet<0)
1001     throw MEDEXCEPTION("MEDMEM::SUPPORT::sub");
1002   else if(sizeRet==0)
1003   {
1004     return 0;
1005   }
1006   if(idsToSuppress==0)
1007   {
1008     ret=new list<int>;
1009     for(int l=0;l<size;l++)
1010       ret->push_back(start+l);
1011     return ret;
1012   }
1013   ret=new list<int>;
1014   int *temp=new int[lgthIdsToSuppress+1];
1015   memcpy(temp,idsToSuppress,sizeof(int)*lgthIdsToSuppress);
1016   temp[lgthIdsToSuppress] = -1;
1017   qsort(temp,lgthIdsToSuppress,sizeof(int),compareId);
1018   int k=0;
1019   for(int i=start;i<=end;i++)
1020     if(temp[k]!=i)
1021       ret->push_back(i);
1022     else
1023       k++;
1024   delete [] temp;
1025   return ret;
1026 }
1027
1028 /*!
1029   performs a common operation : Sub builds a sorted int array that is obtained by supression of all ids contained
1030   in array defined by (idsToSuppress,lgthIdsToSuppress) from array [start ... end]
1031   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 !
1032 */
1033 list<int> *MEDMEM::SUPPORT::sub(const int *ids,int lgthIds,const int *idsToSuppress,int lgthIdsToSuppress)
1034 {
1035   list<int> *ret;
1036   int i,j=0;
1037   if(lgthIds<0)
1038     throw MEDEXCEPTION("MEDMEM::SUPPORT::sub");
1039   else if(lgthIds==0)
1040     return 0;
1041   ret=new list<int>;
1042   int *temp1=new int[lgthIds];
1043   memcpy(temp1,ids,sizeof(int)*lgthIds);
1044   qsort(temp1,lgthIds,sizeof(int),compareId);
1045   int *temp2=new int[lgthIdsToSuppress];
1046   memcpy(temp2,idsToSuppress,sizeof(int)*lgthIdsToSuppress);
1047   qsort(temp2,lgthIdsToSuppress,sizeof(int),compareId);
1048   for(i=0;i<lgthIds;)
1049   {
1050     if(j>=lgthIdsToSuppress)
1051       ret->push_back(temp1[i++]);
1052     else if(temp1[i]>temp2[j])
1053       j++;
1054     else if(temp1[i]<temp2[j])
1055       ret->push_back(temp1[i++]);
1056     else
1057       i++;
1058   }
1059   delete [] temp1;
1060   delete [] temp2;
1061   return ret;
1062 }
1063
1064 /*!
1065   returns a new SUPPORT (responsability to caller to destroy it)
1066   that is the complement to "this", lying on the same entity than "this".
1067 */
1068 SUPPORT *MEDMEM::SUPPORT::getComplement() const
1069 {
1070   SUPPORT *ret;
1071   const int nbOfElt=_mesh->getNumberOfElements(_entity,MED_ALL_ELEMENTS);
1072   int nbOfEltInSupp=getNumberOfElements(MED_ALL_ELEMENTS);
1073   if(_isOnAllElts || nbOfElt==nbOfEltInSupp)
1074   {
1075     ret=new SUPPORT;
1076     ret->setMesh(_mesh);
1077     ret->setEntity(_entity);
1078     string name="Complement of ";
1079     name+=_name;
1080     ret->setName(name);
1081     return ret;
1082   }
1083   const int *nbs=_number->getValue();
1084   list<int> *ids=sub(1,nbOfElt,nbs,nbOfEltInSupp);
1085   if(_entity==MED_NODE)
1086     ret=_mesh->buildSupportOnNodeFromElementList(*ids,_entity);
1087   else
1088     ret=_mesh->buildSupportOnElementsFromElementList(*ids,_entity);
1089   delete ids;
1090   return ret;
1091 }
1092
1093 /*!
1094   returns a new support the user should delete.
1095 */
1096 SUPPORT *MEDMEM::SUPPORT::substract(const SUPPORT& other) const throw (MEDEXCEPTION)
1097 {
1098   const char * LOC = "SUPPORT *MEDMEM::subtract(const SUPPORT * other) : ";
1099   BEGIN_OF_MED(LOC);
1100   SUPPORT *ret;
1101   if (_entity!=other.getEntity())
1102     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Entities are different !"));
1103   if(!(*_mesh == *other.getMesh()))
1104     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Mesh are different !"));
1105   if(other._isOnAllElts)
1106   {
1107     ret=new SUPPORT;
1108     ret->setMesh(_mesh);
1109     ret->setEntity(_entity);
1110     return ret;
1111   }
1112   if(_isOnAllElts)
1113     return other.getComplement();
1114   int nbOfEltInThis=getNumberOfElements(MED_ALL_ELEMENTS);
1115   const int *nbsThis=_number->getValue();
1116   int nbOfEltInOther=other.getNumberOfElements(MED_ALL_ELEMENTS);
1117   const int *nbsOther=other._number->getValue();
1118   list<int> *ids=sub(nbsThis,nbOfEltInThis,nbsOther,nbOfEltInOther);
1119   if(_entity==MED_NODE)
1120     ret=_mesh->buildSupportOnNodeFromElementList(*ids,_entity);
1121   else
1122     ret=_mesh->buildSupportOnElementsFromElementList(*ids,_entity);
1123   delete ids;
1124   return ret;
1125   END_OF_MED(LOC);
1126 }
1127
1128 /*!
1129   returns a new support the user has to delete. Entity is either MED_NODE to obtain node elements lying on boundary of "this"
1130   or MED_FACE,MED_EDGE (depends on the this->_mesh dimension).
1131 */
1132 SUPPORT *MEDMEM::SUPPORT::getBoundaryElements(medEntityMesh Entity) const throw (MEDEXCEPTION)
1133 {
1134   const char * LOC = "SUPPORT *MEDMEM::SUPPORT::getBoundaryElements(medEntityMesh Entity) : ";
1135   BEGIN_OF_MED(LOC);
1136   int spaceDimension=_mesh->getSpaceDimension();
1137   medEntityMesh baseEntity=Entity;
1138   if (spaceDimension == 3)
1139     {
1140       if (Entity!=MED_FACE)
1141         {
1142           if(Entity==MED_NODE)
1143             baseEntity=MED_FACE;
1144           else
1145             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<Entity<<" !"));
1146         }
1147     }
1148   if (spaceDimension == 2)
1149     {
1150       if (Entity!=MED_EDGE)
1151         {
1152           if(Entity==MED_NODE)
1153             baseEntity=MED_EDGE;
1154           else
1155             throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<Entity<<" !"));
1156         }
1157     }
1158   if(_isOnAllElts)
1159     return _mesh->getBoundaryElements(Entity);
1160
1161   const MESH* mesh = _mesh->convertInMESH();
1162   const int * myConnectivityValue=mesh->getReverseConnectivity(MED_DESCENDING);
1163   const int * myConnectivityIndex=mesh->getReverseConnectivityIndex(MED_DESCENDING);
1164   int numberOf=mesh->getNumberOfElements(baseEntity,MED_ALL_ELEMENTS);
1165   const int *ids=_number->getValue();
1166   set<int> idsSet(ids,ids+_totalNumberOfElements);
1167   list<int> myElementsList;
1168   for (int i=0;i<numberOf;i++)
1169   {
1170     int nbOfDP1EntitySharing=0;
1171     if(idsSet.find(myConnectivityValue[myConnectivityIndex[i]-1])!=idsSet.end())
1172       nbOfDP1EntitySharing++;
1173     if(idsSet.find(myConnectivityValue[myConnectivityIndex[i]])!=idsSet.end())
1174       nbOfDP1EntitySharing++;
1175     if(nbOfDP1EntitySharing==1)
1176       myElementsList.push_back(i+1);
1177   }
1178   mesh->removeReference();
1179
1180   if(Entity==MED_NODE)
1181   {
1182     return mesh->buildSupportOnNodeFromElementList(myElementsList,baseEntity);
1183   }
1184   else
1185   {
1186     return mesh->buildSupportOnElementsFromElementList(myElementsList,baseEntity);
1187   }
1188 }
1189
1190 /*!
1191  * \brief Builds a nodal SUPPORT basing on nodes of this one
1192  */
1193
1194 SUPPORT* SUPPORT::buildSupportOnNode() const throw (MEDEXCEPTION)
1195 {
1196   const char * LOC = "SUPPORT *MEDMEM::SUPPORT::buildSupportOnNode() : ";
1197   if ( !getMesh() )
1198     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"This SUPPORT has no mesh"));
1199
1200   SUPPORT* nodalSupport = 0;
1201   if ( isOnAllElements() )
1202     {
1203       nodalSupport = const_cast<SUPPORT*>( getMesh()->getSupportOnAll( MED_NODE ));
1204       nodalSupport->addReference();
1205     }
1206   else
1207     {
1208       if ( !_numberOfElements )
1209         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No element numbers in a partial support"));
1210       string name("Support On Node built from ");
1211       name += getName();
1212
1213       nodalSupport = new SUPPORT;
1214       nodalSupport->setMesh( getMesh());
1215       nodalSupport->setName( name );
1216       nodalSupport->setEntity( MED_NODE );
1217       nodalSupport->setEntity( getEntity() );
1218
1219       const int * nums = _number->getValue();
1220       list<int> elems( nums, nums + _totalNumberOfElements );
1221       getMesh()->fillSupportOnNodeFromElementList( elems, nodalSupport );
1222     }
1223   return nodalSupport;
1224 }
1225
1226 /*!
1227   Method that fills this and updates all its attributes in order to lye on the the listOfNode.
1228 */
1229 void MEDMEM::SUPPORT::fillFromNodeList(const list<int>& listOfNode) throw (MEDEXCEPTION)
1230 {
1231   setEntity(MED_NODE);
1232   clearDataOnNumbers();
1233   int size=listOfNode.size();
1234   int totalNbInMesh=_mesh->getNumberOfElements(_entity,MED_ALL_ELEMENTS);
1235   if(totalNbInMesh==size)
1236   {
1237     _isOnAllElts=true;
1238     update();
1239     return;
1240   }
1241   else
1242     _isOnAllElts=false;
1243   int numberOfGeometricType=1;
1244   medGeometryElement* geometricType=new medGeometryElement[1];
1245   geometricType[0]=MED_NONE;
1246   int *numberOfElements=new int[1];
1247   numberOfElements[0]=size;
1248   int *mySkyLineArrayIndex=new int[2];
1249   mySkyLineArrayIndex[0]=1;
1250   mySkyLineArrayIndex[1]=1+numberOfElements[0];
1251   int *tab=new int[numberOfElements[0]];
1252   int i=0;
1253   for(list<int>::const_iterator iter2=listOfNode.begin();iter2!=listOfNode.end();iter2++)
1254     tab[i++]=*iter2;
1255   MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(1,numberOfElements[0],mySkyLineArrayIndex,tab,true);
1256   setNumberOfGeometricType(numberOfGeometricType);
1257   setGeometricType(geometricType);
1258   setNumberOfElements(numberOfElements);
1259   setNumber(mySkyLineArray);
1260
1261   delete[] numberOfElements;
1262   delete[] geometricType;
1263 }
1264
1265 /*
1266   Method created to factorize code. This method fills the current SUPPORT on entity 'entity' containing all the entities contained in
1267   elements 'listOfElt' of entity 'entity'. Warning this method should be called after both the attributes this->_mesh and this->_entity are correctly set.
1268 */
1269 void MEDMEM::SUPPORT::fillFromElementList(const list<int>& listOfElt) throw (MEDEXCEPTION)
1270 {
1271   clearDataOnNumbers();
1272   int size=listOfElt.size();
1273   int totalNbInMesh=_mesh->getNumberOfElements(_entity,MED_ALL_ELEMENTS);
1274   if(totalNbInMesh==size)
1275   {
1276     _isOnAllElts=true;
1277     update();
1278     return;
1279   }
1280   else
1281     _isOnAllElts=false;
1282   // Well, we must know how many geometric type we have found
1283   int * myListArray = new int[size] ;
1284   int id = 0 ;
1285   list<int>::const_iterator myElementsListIt ;
1286   for (myElementsListIt=listOfElt.begin();myElementsListIt!=listOfElt.end();myElementsListIt++)
1287     myListArray[id++]=(*myElementsListIt) ;
1288   int numberOfGeometricType ;
1289   medGeometryElement* geometricType ;
1290   int * numberOfElements ;
1291   int * mySkyLineArrayIndex ;
1292
1293   int numberOfType = _mesh->getNumberOfTypes(_entity) ;
1294   if (numberOfType == 1)
1295   {
1296     numberOfGeometricType = 1 ;
1297     geometricType = new medGeometryElement[1] ;
1298     geometricType[0] = _mesh->getTypes(_entity)[0];
1299     numberOfElements = new int[1] ;
1300     numberOfElements[0] = size ;
1301     mySkyLineArrayIndex = new int[2] ;
1302     mySkyLineArrayIndex[0]=1 ;
1303     mySkyLineArrayIndex[1]=1+size ;
1304   }
1305   else // hemmm
1306   {
1307     map<medGeometryElement,int> theType ;
1308     for (myElementsListIt=listOfElt.begin();myElementsListIt!=listOfElt.end();myElementsListIt++)
1309     {
1310       medGeometryElement myType = _mesh->getElementType(_entity,*myElementsListIt) ;
1311       if (theType.find(myType) != theType.end() )
1312         theType[myType]+=1 ;
1313       else
1314         theType[myType]=1 ;
1315     }
1316     numberOfGeometricType = theType.size() ;
1317     geometricType = new medGeometryElement[numberOfGeometricType] ;
1318     numberOfElements = new int[numberOfGeometricType] ;
1319     mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
1320     int index = 0 ;
1321     mySkyLineArrayIndex[0]=1 ;
1322     map<medGeometryElement,int>::iterator theTypeIt ;
1323     for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++)
1324     {
1325       geometricType[index] = (*theTypeIt).first ;
1326       numberOfElements[index] = (*theTypeIt).second ;
1327       mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfElements[index] ;
1328       index++ ;
1329     }
1330   }
1331   MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray,/*shallowCopy=*/true) ;
1332   setNumberOfGeometricType(numberOfGeometricType) ;
1333   setGeometricType(geometricType) ;
1334   setNumberOfElements(numberOfElements) ;
1335   setNumber(mySkyLineArray) ;
1336
1337   delete[] numberOfElements;
1338   delete[] geometricType;
1339 }
1340
1341 /*! \ifdef SUPPORT_advanced
1342
1343 \brief creates a MESH that contains only the elements in the current support.
1344
1345 The output mesh has no group, nor elements of connectivity lesser than that of the present support. Nodes are renumbered so that they are numberd from 1 to N in the new mesh. The order of the elements in the new mesh corresponds to that of the elements in the original support.
1346 */
1347 MESH* SUPPORT::makeMesh() const
1348 {
1349   const char* LOC = "SUPPORT::makeMesh(): ";
1350   if ( !_mesh )
1351     throw MED_EXCEPTION(STRING(LOC)<<" NULL mesh in support");
1352   if ( _entity == MED_NODE )
1353     throw MED_EXCEPTION(STRING(LOC)<<" unavailable for support on nodes");
1354
1355   //Creating the new mesh
1356
1357   MESHING* newmesh = new MESHING();
1358   newmesh->setName( STRING("MeshFromSupport_") << getName() );
1359
1360   // set types info
1361   const medGeometryElement* types = getTypes();
1362   int nb_types = _numberOfGeometricType;
1363   newmesh->setNumberOfTypes   ( nb_types, MED_CELL );
1364   newmesh->setTypes           ( types, MED_CELL );
1365   newmesh->setNumberOfElements( _numberOfElements, MED_CELL);
1366
1367   // browsing through elements to create a mapping between
1368   // the new nodes and the old nodes and to create nodal connectivity
1369
1370   const MESH* mesh = _mesh->convertInMESH();
1371   const medGeometryElement* all_mesh_types = mesh->getTypes( _entity );
1372   const int *                    num_index = mesh->getGlobalNumberingIndex( _entity );
1373
1374   map<int,int> oldnodes; // map old to new nodes
1375   if ( types[nb_types-1] == MED_POLYHEDRA )
1376     oldnodes.insert( make_pair( -1, -1 )); // for face separators
1377   int newid=1;
1378   for (int itype=0; itype < _numberOfGeometricType;itype++)
1379   {
1380     medGeometryElement type = types[itype];
1381     int nbelems = getNumberOfElements(type);
1382
1383     // get connectivity info
1384     int shift = 1; // to pass from elem number to array index
1385     const int* conn = mesh->getConnectivity(MED_NODAL,_entity,MED_ALL_ELEMENTS);
1386     const int* index = mesh->getConnectivityIndex(MED_NODAL,_entity);
1387     int t = 0;
1388     while ( type != all_mesh_types[t] ) ++t;
1389     shift+= num_index[ t ] - num_index[0];
1390     index+= num_index[ t ] - num_index[0];
1391
1392     // make and set new connectivity
1393     if ( _isOnAllElts  && _entity == MED_CELL )
1394     {
1395       newmesh->setConnectivity( MED_CELL, type, conn, index );
1396     }
1397     else // partial support or support of sub-entities
1398     {
1399       vector<int> new_conn, new_index;
1400       new_conn.reserve ( nbelems * (type % 100) );
1401       new_index.reserve( nbelems + 1 );
1402       new_index.push_back(1);
1403
1404       const int * nums = _isOnAllElts ? 0 : getNumber( type );
1405
1406       for (int i=0; i<nbelems; i++)
1407       {
1408         int ielem = nums ? nums[i]-shift : i;
1409         const int* elem_node = conn + index[ ielem   ] - 1;
1410         const int* nodes_end = conn + index[ ielem+1 ] - 1;
1411         for ( ; elem_node < nodes_end; ++elem_node )
1412         {
1413           // make new connectivity
1414           map<int,int>::iterator old_new = oldnodes.insert(make_pair( *elem_node, newid )).first;
1415           new_conn.push_back( old_new->second );
1416           if ( old_new->second == newid )
1417             newid++;
1418         }
1419         new_index.push_back( new_index.back() + index[ ielem+1 ] - index[ ielem ] );
1420       }
1421       // set new connectivity
1422       newmesh->setConnectivity( MED_CELL, type, & new_conn[0], & new_index[0] );
1423     }
1424   }
1425
1426   //definition of coordinates
1427
1428   int nb_nodes, spacedim = mesh->getSpaceDimension();
1429   const double*oldcoords = mesh->getCoordinates(MED_FULL_INTERLACE);
1430   PointerOf<double> newcoords;
1431
1432   if ( _isOnAllElts && _entity == MED_CELL )
1433   {
1434     nb_nodes = mesh->getNumberOfNodes();
1435     newcoords.set( oldcoords );
1436   }
1437   else
1438   {
1439     nb_nodes = oldnodes.size();
1440     newcoords.set( nb_nodes * spacedim);
1441     for (std::map<int,int>::const_iterator iter=oldnodes.begin(); iter!=oldnodes.end();iter++)
1442       std::copy( oldcoords+(iter->first-1 )*spacedim,
1443                  oldcoords+(iter->first   )*spacedim,
1444                  newcoords+(iter->second-1)*spacedim);
1445   }
1446   newmesh->setCoordinates(spacedim, nb_nodes, newcoords,
1447                           mesh->getCoordinatesSystem(), MED_FULL_INTERLACE);
1448   newmesh->setCoordinatesNames ( mesh->getCoordinatesNames() );
1449   newmesh->setCoordinatesUnits ( mesh->getCoordinatesUnits() );
1450
1451   mesh->removeReference();
1452
1453   return newmesh;
1454 }
1455 /*!
1456   @}
1457 */
1458
1459 /*! set the reference _mesh to Mesh */
1460 //--------------------------------------
1461 void SUPPORT::setMesh(const GMESH *Mesh) const
1462   //--------------------------------------
1463 {
1464   if(_mesh!=Mesh)
1465     {
1466       if(_mesh)
1467         _mesh->removeReference();
1468       _mesh=Mesh;
1469       _meshName = "";
1470       if(_mesh)
1471         _mesh->addReference();
1472     }
1473 }
1474
1475 /*! returns the mesh name  */
1476 //------------------------------------
1477 string SUPPORT::getMeshName() const
1478   //------------------------------------
1479 {
1480   if (_mesh)
1481     return _mesh->getName();
1482   else
1483     return _meshName;
1484 }
1485
1486 /*!\if MEDMEM_ug
1487 \addtogroup SUPPORT_query
1488 @{
1489 \endif
1490 */
1491
1492 /*!
1493   This method returns the number of all elements of the type GeometricType.
1494
1495   If isOnAllElements is false, it returns the number of elements in the
1496   support otherwise it returns number of elements in the mesh.
1497
1498   Example : number of MED_TRIA3 or MED_ALL_ELEMENTS elements
1499   in support.
1500
1501   Note : If SUPPORT is defined on MED_NODE, use MED_ALL_ELEMENTS as
1502          medGeometryElement GeometricType and it will return the number
1503          of nodes in the support (or in the mesh).
1504 */
1505 //-----------------------------------------------------------------------------
1506 int SUPPORT::getNumberOfElements(MED_EN::medGeometryElement GeometricType) const
1507   throw (MEDEXCEPTION)
1508 //-----------------------------------------------------------------------------
1509 {
1510   if (GeometricType==MED_EN::MED_ALL_ELEMENTS)
1511     return _totalNumberOfElements;
1512   for (int i=0;i<_numberOfGeometricType;i++)
1513     if (_geometricType[i]==GeometricType)
1514       return ( _totalNumberOfElements < 1 ) ? 0 : _numberOfElements[i];
1515   throw MEDEXCEPTION("Support::getNumberOfElements : Geometric type not found !") ;
1516 }
1517
1518   /*! Returns the total number of elements in the support. */
1519 //-----------------------------------------------------------------------------
1520 const int * SUPPORT::getNumberOfElements() const throw (MEDEXCEPTION) {
1521 //-----------------------------------------------------------------------------
1522   return _numberOfElements;
1523 }
1524
1525 /*!
1526   Returns index of element number.
1527
1528   Note : See getConnectivityIndex for details.
1529 */
1530 //-------------------------------------------
1531 const int * SUPPORT::getNumberIndex() const
1532 //-------------------------------------------
1533   throw (MEDEXCEPTION)
1534 {
1535   /*  issue 0021167: [CEA 448] Supports management on all elements
1536   if (_isOnAllElts)
1537     throw MEDEXCEPTION("Support::getNumberIndex : Not defined, support is on all entity !") ;
1538   */
1539   if ( !_number )
1540     throw MEDEXCEPTION("Support::getNumberIndex : numbers not set !") ;
1541   return _number->getIndex() ;
1542 }
1543 /*! \if MEDMEM_ug
1544 @}
1545 \endif */
1546
1547 //---------------------------------------------------------------------
1548 MEDSKYLINEARRAY * SUPPORT::getnumber() const
1549   throw (MEDEXCEPTION)
1550 //---------------------------------------------------------------------
1551 {
1552   if (_number==NULL)
1553     throw MEDEXCEPTION("Support::getnumber : Not defined !") ;
1554   return _number ;
1555 }
1556
1557 //---------------------------------------------------------------------
1558 MEDSKYLINEARRAY * SUPPORT::getnumberFromFile() const
1559   throw (MEDEXCEPTION)
1560 //---------------------------------------------------------------------
1561 {
1562   if (_number_fromfile==NULL)
1563     throw MEDEXCEPTION("Support::getnumberFromFile : Not defined !") ;
1564   return _number_fromfile ;
1565 }
1566
1567 /*!
1568   Returns an array which contains all number of given medGeometryElement.
1569
1570   Numbering is global, ie numbers are bounded by 1 and
1571   GMESH::getNumberOfElement(entity,MED_ALL_ELEMENTS) and not by 1 and
1572   GMESH::getNumberOfElement(entity,geomElement).
1573
1574   Note : If SUPPORT is defined on MED_NODE, use MED_NONE
1575   medGeometryElement type.
1576 */
1577 //---------------------------------------------------------------------
1578 const int * SUPPORT::getNumber(MED_EN::medGeometryElement GeometricType) const
1579   throw (MEDEXCEPTION)
1580 //---------------------------------------------------------------------
1581 {
1582   /*  issue 0021167: [CEA 448] Supports management on all elements
1583   if (_isOnAllElts)
1584     throw MEDEXCEPTION("Support::getNumber : Not defined, support is on all entity !") ;
1585   */
1586   if (!_number)
1587   {
1588     if ( _isOnAllElts )
1589       //update();
1590       throw MEDEXCEPTION("Support::getNumber : not updated (update()) SUPPORT on all elements !") ;
1591     else if ( _totalNumberOfElements > 0 )
1592       throw MEDEXCEPTION("Support::getNumber : wrong support, _number not defined !") ;
1593     else
1594       return NULL;
1595   }
1596   if (GeometricType==MED_EN::MED_ALL_ELEMENTS)
1597     return _number->getValue() ;
1598   for (int i=0;i<_numberOfGeometricType;i++)
1599     if (_geometricType[i]==GeometricType)
1600       return _number->getI(i+1) ;
1601   throw MEDEXCEPTION("Support::getNumber : GeometricType not found !") ;
1602 }
1603
1604 //---------------------------------------------------------------------
1605 const int * SUPPORT::getNumberFromFile(MED_EN::medGeometryElement GeometricType) const
1606   throw (MEDEXCEPTION)
1607 //---------------------------------------------------------------------
1608 {
1609 //   if (_isOnAllElts)
1610 //     throw MEDEXCEPTION("Support::getNumberFromFile : Not defined, support is on all entity !") ;
1611   if (GeometricType==MED_EN::MED_ALL_ELEMENTS)
1612     return _number_fromfile->getValue() ;
1613   for (int i=0;i<_numberOfGeometricType;i++)
1614     if (_geometricType[i]==GeometricType)
1615       return _number_fromfile->getI(i+1) ;
1616   throw MEDEXCEPTION("Support::getNumberFromFile : GeometricType not found !") ;
1617 }
1618
1619 /*! set the meshName if there is ni reference _mesh to Mesh */
1620 //--------------------------------------
1621 void SUPPORT::setMeshName(const std::string & meshName)
1622 //--------------------------------------
1623 {
1624   if (_mesh)
1625     throw MEDEXCEPTION("SUPPORT::setMeshName(const string & meshName) : Setting meshName is not possible when an associated mesh is set !") ;
1626
1627   _meshName=meshName;
1628 }
1629
1630 /*! set the attribute _entity to Entity */
1631 //------------------------------------------
1632 void SUPPORT::setEntity(MED_EN::medEntityMesh Entity)
1633 {
1634   _entity=Entity;
1635   // 0021199: [CEA 458] MEDMEM::SUPPORT : geometric type when support is on node
1636   // set geometric type -> MED_NONE
1637   if ( _entity == MED_NODE )
1638     {
1639       _numberOfGeometricType = 1;
1640       const MED_EN::medGeometryElement nodeType = MED_EN::MED_NONE;
1641       _geometricType.set(0);
1642       setGeometricType( &nodeType );
1643     }
1644 }
1645
1646 /*! set the attribute _numberOfGeometricType to NumberOfGeometricType */
1647 //---------------------------------------------------------------------
1648 void SUPPORT::setNumberOfGeometricType(int NumberOfGeometricType)
1649 //---------------------------------------------------------------------
1650 {
1651   _numberOfGeometricType=NumberOfGeometricType;
1652
1653   _geometricType.set(0);
1654   _numberOfElements.set(0);
1655   _profilNames.resize( NumberOfGeometricType, "" );
1656 }
1657
1658 /*! set the attribute _geometricType to geometricType */
1659 //---------------------------------------------------------------------
1660 void SUPPORT::setGeometricType(const MED_EN::medGeometryElement *GeometricType)
1661 //---------------------------------------------------------------------
1662 {
1663   if (!_geometricType)
1664     {
1665       _geometricType.set(_numberOfGeometricType, GeometricType);
1666       // 0021199: [CEA 458] MEDMEM::SUPPORT : geometric type when support is on node
1667       // geometric type must be MED_NONE
1668       if ( _entity == MED_NODE && _numberOfGeometricType == 1 && _geometricType[0] != MED_NONE )
1669         throw MEDEXCEPTION("SUPPORT::setGeometricType(), valid type for MED_NODE is MED_NONE ");
1670     }
1671   if (_profilNames.empty() || _profilNames[0].empty())
1672     {
1673       // giving a default value to profile names
1674       vector<string> prof_names( _numberOfGeometricType);
1675       string name = healName( _name );
1676       for (int itype=0; itype < _numberOfGeometricType; itype++)
1677         {
1678           ostringstream typestr;
1679           typestr<<name<<"_type"<<_geometricType[itype];
1680           prof_names[itype]=typestr.str();
1681         }
1682       _profilNames=prof_names;
1683     }
1684 }
1685
1686 /*!
1687   Set the attribute _numberOfElements to NumberOfElements and
1688   calculate the total number of elements.
1689 */
1690 //----------------------------------------------------------
1691 void SUPPORT::setNumberOfElements(const int *NumberOfElements)
1692 //----------------------------------------------------------
1693 {
1694   if (_numberOfElements == NULL)
1695     {
1696       if (_numberOfGeometricType)
1697         _numberOfElements.set(_numberOfGeometricType,NumberOfElements);
1698       else
1699         _numberOfElements.set(0);
1700     }
1701   _totalNumberOfElements = 0 ;
1702   for (int i=0;i<_numberOfGeometricType;i++)
1703     _totalNumberOfElements+=_numberOfElements[i];
1704 }
1705
1706 /*! set the attribute _number to Number */
1707 //---------------------------------------------------
1708 void SUPPORT::setNumber(MEDSKYLINEARRAY * Number)
1709 //---------------------------------------------------
1710 {
1711   /*  issue 0021167: [CEA 448] Supports management on all elements
1712   if ( _isOnAllElts )
1713     throw MEDEXCEPTION("SUPPORT::setNumber(MEDSKYLINEARRAY * Number) Support is on all elements") ;
1714   */
1715   if (_number != NULL) delete _number ;
1716   _number=Number;
1717 }
1718
1719 /*! set the attribute _number with index and value arrays */
1720 //---------------------------------------------------
1721 void SUPPORT::setNumber(const int * index, const int* value, bool shallowCopy)
1722 //---------------------------------------------------
1723 {
1724   if (_number != NULL) delete _number ;
1725   _number= new MEDSKYLINEARRAY(_numberOfGeometricType,_totalNumberOfElements,index,value,shallowCopy);
1726 }