src/MEDMEM/MEDMEM_Mesh.cxx

   1 //  Copyright (C) 2007-2008  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  File Mesh.cxx
  24 */
  25
  26 #include <math.h>
  27 #include <map>
  28 #include <sstream>
  29
  30 #include "MEDMEM_DriversDef.hxx"
  31 #include "MEDMEM_Field.hxx"
  32 #include "MEDMEM_Mesh.hxx"
  33
  34 #include "MEDMEM_Support.hxx"
  35 #include "MEDMEM_Family.hxx"
  36 #include "MEDMEM_Group.hxx"
  37 #include "MEDMEM_Coordinate.hxx"
  38 #include "MEDMEM_Connectivity.hxx"
  39 #include "MEDMEM_CellModel.hxx"
  40 #include "VolSurfFormulae.hxx"
  41 #include "MEDMEM_InterpolationHighLevelObjects.hxx"
  42 #include "MEDMEM_DriverFactory.hxx"
  43
  44 using namespace std;
  45 using namespace MEDMEM;
  46 using namespace MED_EN;
  47
  48 #define MED_NOPDT -1
  49 #define MED_NONOR -1
  50
  51 // Block defining groups for the MEDMEM_ug documentation
  52 /*!
  53 \defgroup MESH_constructors MESH Constructors
  54
  55 The MESH class provides only two constructors : a copy constructor and
  56 a constructor enabling creation from file reading. The creation of
  57 a user-defined mesh implies the use of the MESHING class.
  58
  59 \defgroup MESH_advanced MESH Advanced features
  60 These functions provide access to high-level manipulation of the meshes, giving
  61 information about the cells or extracting supports from the mesh.
  62
  63 \defgroup MESH_connectivity MESH Connectivity information
  64 These methods are related to the extraction of connectivity information
  65 from the mesh.
  66
  67 \defgroup MESH_nodes MESH Nodes information
  68 These methods are related to the extraction of information about the mesh nodes.
  69
  70 \defgroup MESH_general MESH General information
  71
  72 These methods are related to the retrieval of general information about the mesh.
  73
  74 \defgroup MESH_poly MESH Polygons and Polyhedra information
  75
  76 These methods are specific methods used for retrieving connectivity
  77 information for MED_POLYGON and MED_POLYHEDRON elements.
  78
  79
  80 \defgroup MESH_families Families and Groups handling
  81
  82 The methods described in this section enable the manipulation of families and groups. These
  83 notions define subsets of MED elements in a mesh. They differ because families are non
  84 overlapping (a mesh element is associated to zero or one family)  while groups are more general.
  85
  86 \defgroup MESH_io Mesh I/O
  87 These methods describe how to read and write meshes. Generally speaking, meshes should be read
  88 via a constructor and should be written with the addDriver()/write() methods.
  89
  90 */
  91
  92 // ------- Drivers Management Part
  93
  94 /*! Add a %MESH driver of type %driverTypes (MED_DRIVER, ....) associated with file fileName. The meshname used in the file
  95     is  driverName. addDriver returns an integer handler. */
  96 int MESH::addDriver(driverTypes driverType,
  97                     const string & fileName/*="Default File Name.med"*/,
  98                     const string & driverName/*="Default Mesh Name"*/,
  99                     MED_EN::med_mode_acces access)
 100 {
 101   const char* LOC = "MESH::addDriver(driverTypes driverType, const string & fileName=\"Default File Name.med\",const string & driverName=\"Default Mesh Name\",MED_EN::med_mode_acces access) : ";
 102   BEGIN_OF_MED(LOC);
 103
 104   GENDRIVER * driver;
 105
 106   SCRUTE_MED(driverType);
 107
 108   driver = DRIVERFACTORY::buildDriverForMesh(driverType,fileName,this,
 109                                              driverName,access) ;
 110
 111   _drivers.push_back(driver);
 112
 113   int current = _drivers.size()-1;
 114
 115   _drivers[current]->setMeshName(driverName);
 116
 117   END_OF_MED(LOC);
 118
 119   return current;
 120 }
 121
 122 /*! Add an existing MESH driver. */
 123 int  MESH::addDriver(GENDRIVER & driver)
 124 {
 125   const char* LOC = "MESH::addDriver(GENDRIVER &) : ";
 126   BEGIN_OF_MED(LOC);
 127
 128   // A faire : Vérifier que le driver est de type MESH.
 129
 130   // For the case where driver does not know about me, i.e. has been created through
 131   // constructor witout parameters: create newDriver knowing me and get missing data
 132   // from driver using merge()
 133   //GENDRIVER * newDriver = driver.copy() ;
 134   GENDRIVER* newDriver = DRIVERFACTORY::buildDriverForMesh(driver.getDriverType(),
 135                                                            driver.getFileName(), this,
 136                                                            driver.getMeshName(),
 137                                                            driver.getAccessMode());
 138   _drivers.push_back(newDriver);
 139
 140   int current = _drivers.size()-1;
 141   driver.setId(current);
 142
 143   newDriver->merge( driver );
 144   newDriver->setId( current );
 145
 146   return current;
 147
 148   END_OF_MED(LOC);
 149 }
 150
 151 /*! Remove an existing MESH driver. */
 152 void MESH::rmDriver (int index/*=0*/) {
 153   const char * LOC = "MESH::rmDriver (int index=0): ";
 154   BEGIN_OF_MED(LOC);
 155
 156   if (index >= 0 && index < _drivers.size() && _drivers[index]) {
 157     delete _drivers[index];
 158     _drivers[index] = 0;
 159      MESSAGE_MED ("detruire");
 160   }
 161   else
 162     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
 163                                      << "The index given is invalid, index must be between  0 and  |"
 164                                      << _drivers.size()
 165                                      )
 166                           );
 167
 168   END_OF_MED(LOC);
 169
 170 };
 171
 172 // ------ End of Drivers Management Part
 173
 174
 175 void MESH::init()
 176 {
 177
 178   const char* LOC = "MESH::init(): ";
 179   BEGIN_OF_MED(LOC);
 180
 181   _name = "NOT DEFINED"; // A POSITIONNER EN FCT DES IOS ?
 182
 183   _coordinate   = (COORDINATE   *) NULL;
 184   _connectivity = (CONNECTIVITY *) NULL;
 185
 186   _spaceDimension        =          MED_INVALID; // 0 ?!?
 187   _meshDimension         =          MED_INVALID;
 188   _numberOfNodes         =          MED_INVALID;
 189
 190   _isAGrid = false;
 191
 192   _arePresentOptionnalNodesNumbers = 0;
 193
 194   END_OF_MED(LOC);
 195 };
 196
 197 /*! Create an empty MESH. */
 198 MESH::MESH():_coordinate(NULL),_connectivity(NULL), _isAGrid(false) {
 199   init();
 200 };
 201
 202 /*! \if MEDMEM_ug
 203   \addtogroup MESH_constructors
 204 @{
 205 \endif
 206 */
 207 /*!
 208 Copy constructor
 209 */
 210 MESH::MESH(MESH &m)
 211 {
 212   _name=m._name;
 213   _isAGrid = m._isAGrid;
 214
 215   if (m._coordinate != NULL)
 216     _coordinate = new COORDINATE(* m._coordinate);
 217   else
 218     _coordinate = (COORDINATE *) NULL;
 219
 220   if (m._connectivity != NULL)
 221     _connectivity = new CONNECTIVITY(* m._connectivity);
 222   else
 223     _connectivity = (CONNECTIVITY *) NULL;
 224
 225   _spaceDimension = m._spaceDimension;
 226   _meshDimension = m._meshDimension;
 227   _numberOfNodes = m._numberOfNodes;
 228
 229   _familyNode = m._familyNode;
 230   for (int i=0; i<(int)m._familyNode.size(); i++)
 231     {
 232       _familyNode[i] = new FAMILY(* m._familyNode[i]);
 233       _familyNode[i]->setMeshDirectly(this);
 234     }
 235
 236   _familyCell = m._familyCell;
 237   for (int i=0; i<(int)m._familyCell.size(); i++)
 238     {
 239       _familyCell[i] = new FAMILY(* m._familyCell[i]);
 240       _familyCell[i]->setMeshDirectly(this);
 241     }
 242
 243   _familyFace = m._familyFace;
 244   for (int i=0; i<(int)m._familyFace.size(); i++)
 245     {
 246       _familyFace[i] = new FAMILY(* m._familyFace[i]);
 247       _familyFace[i]->setMeshDirectly(this);
 248     }
 249
 250   _familyEdge = m._familyEdge;
 251   for (int i=0; i<(int)m._familyEdge.size(); i++)
 252     {
 253       _familyEdge[i] = new FAMILY(* m._familyEdge[i]);
 254       _familyEdge[i]->setMeshDirectly(this);
 255     }
 256
 257   _groupNode = m._groupNode;
 258   for (int i=0; i<(int)m._groupNode.size(); i++)
 259     {
 260       _groupNode[i] = new GROUP(* m._groupNode[i]);
 261       _groupNode[i]->setMeshDirectly(this);
 262     }
 263
 264   _groupCell = m._groupCell;
 265   for (int i=0; i<(int)m._groupCell.size(); i++)
 266     {
 267       _groupCell[i] = new GROUP(* m._groupCell[i]);
 268       _groupCell[i]->setMeshDirectly(this);
 269     }
 270
 271   _groupFace = m._groupFace;
 272   for (int i=0; i<(int)m._groupFace.size(); i++)
 273     {
 274       _groupFace[i] = new GROUP(* m._groupFace[i]);
 275       _groupFace[i]->setMeshDirectly(this);
 276     }
 277
 278   _groupEdge = m._groupEdge;
 279   for (int i=0; i<(int)m._groupEdge.size(); i++)
 280     {
 281       _groupEdge[i] = new GROUP(* m._groupEdge[i]);
 282       _groupEdge[i]->setMeshDirectly(this);
 283     }
 284
 285   //_drivers = m._drivers;  //Recopie des drivers?
 286 }
 287
 288 /*!
 289 \if MEDMEM_ug
 290 @}
 291 \endif
 292 */
 293
 294 MESH::~MESH() {
 295
 296   MESSAGE_MED("MESH::~MESH() : Destroying the Mesh");
 297   if (_coordinate != ((COORDINATE *) NULL)) delete _coordinate ;
 298   if (_connectivity != ((CONNECTIVITY *) NULL)) delete _connectivity ;
 299   int size ;
 300   size = _familyNode.size() ;
 301   for (int i=0;i<size;i++)
 302     delete _familyNode[i] ;
 303   size = _familyCell.size() ;
 304   for (int i=0;i<size;i++)
 305     delete _familyCell[i] ;
 306   size = _familyFace.size() ;
 307   for (int i=0;i<size;i++)
 308     delete _familyFace[i] ;
 309   size = _familyEdge.size() ;
 310   for (int i=0;i<size;i++)
 311     delete _familyEdge[i] ;
 312   size = _groupNode.size() ;
 313   for (int i=0;i<size;i++)
 314     delete _groupNode[i] ;
 315   size = _groupCell.size() ;
 316   for (int i=0;i<size;i++)
 317     delete _groupCell[i] ;
 318   size = _groupFace.size() ;
 319   for (int i=0;i<size;i++)
 320     delete _groupFace[i] ;
 321   size = _groupEdge.size() ;
 322   for (int i=0;i<size;i++)
 323     delete _groupEdge[i] ;
 324
 325   map<medEntityMesh,SUPPORT*>::iterator it = _entitySupport.begin();
 326   for(;it!=_entitySupport.end();it++)
 327     if((*it).second != NULL)
 328       delete (*it).second;
 329
 330   MESSAGE_MED("In this object MESH there is(are) " << _drivers.size() << " driver(s)");
 331
 332   for (unsigned int index=0; index < _drivers.size(); index++ )
 333     {
 334       SCRUTE_MED(_drivers[index]);
 335       if ( _drivers[index] != NULL) delete _drivers[index];
 336     }
 337
 338 }
 339
 340
 341 /*! \if MEDMEM_ug
 342 \addtogroup MESH_poly
 343 @{
 344 \endif
 345 */
 346 /*!
 347   Method equivalent to getNumberOfTypes except that it includes not only classical Types but polygons/polyhedra also.
 348  */
 349 int MESH::getNumberOfTypesWithPoly(MED_EN::medEntityMesh Entity) const
 350 {
 351   if(_connectivity!= NULL)
 352     return _connectivity->getNumberOfTypesWithPoly(Entity);
 353   throw MEDEXCEPTION(LOCALIZED("MESH::getNumberOfTypesWithPoly( medEntityMesh ) : Connectivity not defined !"));
 354 }
 355
 356 /*
 357   Method equivalent to getTypesWithPoly except that it includes not only classical Types but polygons/polyhedra also.
 358   WARNING the returned array MUST be deallocated.
 359  */
 360 MED_EN::medGeometryElement * MESH::getTypesWithPoly(MED_EN::medEntityMesh Entity) const
 361 {
 362   if (Entity == MED_EN::MED_NODE)
 363     throw MEDEXCEPTION(LOCALIZED("MESH::getTypes( medEntityMesh ) : No medGeometryElement with MED_NODE entity !"));
 364   if (_connectivity != NULL)
 365     return _connectivity->getGeometricTypesWithPoly(Entity);
 366   throw MEDEXCEPTION(LOCALIZED("MESH::getTypes( medEntityMesh ) : Connectivity not defined !"));
 367 }
 368
 369 /*
 370   Method equivalent to getNumberOfElementsWithPoly except that it includes not only classical Types but polygons/polyhedra also.
 371  */
 372 int MESH::getNumberOfElementsWithPoly(MED_EN::medEntityMesh Entity, MED_EN::medGeometryElement Type) const
 373 {
 374   if(Type==MED_POLYGON || Type==MED_POLYHEDRA)
 375     {
 376       int nbOfPolygs=_connectivity->getNumberOfElementOfPolyType(Entity);
 377       return nbOfPolygs;
 378     }
 379   else if(Type==MED_ALL_ELEMENTS)
 380     {
 381       int nbOfClassicalTypes=getNumberOfElements(Entity,MED_ALL_ELEMENTS);
 382       int nbOfClassicalTypes2=_connectivity->getNumberOfElementOfPolyType(Entity);
 383       return nbOfClassicalTypes+nbOfClassicalTypes2;
 384     }
 385   else
 386     return getNumberOfElements(Entity,Type);
 387 }
 388
 389 /*! \if MEDMEM_ug
 390 @}
 391 \endif*/
 392 bool MESH::existConnectivityWithPoly(MED_EN::medConnectivity ConnectivityType,
 393                                      MED_EN::medEntityMesh Entity) const
 394 {
 395   if (_connectivity==(CONNECTIVITY*)NULL)
 396     throw MEDEXCEPTION("MESH::existConnectivity(medConnectivity,medEntityMesh) : no connectivity defined !");
 397   return _connectivity->existConnectivityWithPoly(ConnectivityType,Entity);
 398 }
 399
 400 MESH & MESH::operator=(const MESH &m)
 401 {
 402   const char* LOC = "MESH & MESH::operator=(const MESH &m) : ";
 403   BEGIN_OF_MED(LOC);
 404
 405   MESSAGE_MED(PREFIX_MED <<"Not yet implemented, operating on the object " << m);
 406   //  A FAIRE.........
 407
 408   // ATTENTION CET OPERATEUR DE RECOPIE EST DANGEREUX POUR LES
 409   // POINTEURS : ex : nodal_connectivity ???? EXPRES ????
 410 //    _drivers = m._drivers;
 411
 412 //    space_dimension=m.space_dimension;
 413 //    mesh_dimension=m.mesh_dimension;
 414
 415 //    nodes_count=m.nodes_count;
 416 //    coordinates=m.coordinates;
 417 //    coordinates_name=m.coordinates_name ;
 418 //    coordinates_unit=m.coordinates_unit ;
 419 //    nodes_numbers=m.nodes_numbers ;
 420
 421 //    cells_types_count=m.cells_types_count;
 422 //    cells_type=m.cells_type;
 423 //    cells_count=m.cells_count;
 424 //    nodal_connectivity=m.nodal_connectivity;
 425
 426 //    nodes_families_count=m.nodes_families_count;
 427 //    nodes_Families=m.nodes_Families;
 428
 429 //    cells_families_count=m.cells_families_count;
 430 //    cells_Families=m.cells_Families;
 431
 432 //    maximum_cell_number_by_node = m.maximum_cell_number_by_node;
 433 //    if (maximum_cell_number_by_node > 0)
 434 //      {
 435 //        reverse_nodal_connectivity = m.reverse_nodal_connectivity;
 436 //        reverse_nodal_connectivity_index = m.reverse_nodal_connectivity_index ;
 437 //      }
 438   END_OF_MED(LOC);
 439
 440   return *this;
 441 }
 442
 443 bool MESH::operator==(const MESH& other) const
 444 {
 445   const char* LOC = "MESH::operator==";
 446   BEGIN_OF_MED(LOC);
 447   return this==&other;
 448 }
 449
 450 /*!\if MEDMEM_ug
 451 \addtogroup MESH_constructors
 452 @{
 453 \endif
 454 */
 455 /*!
 456  Creates a %MESH object using a %MESH driver of type %driverTypes (MED_DRIVER, GIBI_DRIVER, ...) associated with file \a fileName. As several meshes can coexist in the same file (notably in MED files) , the constructor takes a third argument that specifies the name of the mesh.
 457 The constructor will throw an exception if the file does not exist, has no reading permissions or if the mesh does not exist in the file.
 458 */
 459 MESH::MESH(driverTypes driverType, const string &  fileName/*=""*/, const string &  driverName/*=""*/) throw (MEDEXCEPTION)
 460 {
 461   const char * LOC = "MESH::MESH(driverTypes driverType, const string &  fileName="", const string &  driverName/="") : ";
 462   BEGIN_OF_MED(LOC);
 463
 464   int current;
 465
 466   init();
 467   GENDRIVER *myDriver=DRIVERFACTORY::buildDriverForMesh(driverType,fileName,this,driverName,RDONLY);
 468   current = addDriver(*myDriver);
 469   delete myDriver;
 470   _drivers[current]->open();
 471   _drivers[current]->read();
 472   _drivers[current]->close();
 473
 474   END_OF_MED(LOC);
 475 };
 476 /*!\if MEDMEM_ug
 477   @}
 478 \endif
 479 */
 480
 481 /*!
 482 \addtogroup MESH_general
 483 @{
 484 */
 485 /*!
 486 Returns true if mesh \a other has same
 487 coordinates (to 1E-15 precision ) and same connectivity as the calling object.
 488 Information like name or description is not taken into account
 489 for the comparison.
 490 */
 491
 492 bool MESH::deepCompare(const MESH& other) const
 493 {
 494   int size1=getSpaceDimension()*getNumberOfNodes();
 495   int size2=other.getSpaceDimension()*other.getNumberOfNodes();
 496   if(size1!=size2)
 497     return false;
 498
 499   const COORDINATE* CRD = other.getCoordinateptr();
 500   if( (!CRD && _coordinate) || (CRD && !(_coordinate)) ) {
 501     return false;
 502   }
 503
 504   bool ret=true;
 505   if( _coordinate ) {
 506     const double* coord1=getCoordinates(MED_FULL_INTERLACE);
 507     const double* coord2=other.getCoordinates(MED_FULL_INTERLACE);
 508     for(int i=0;i<size1 && ret;i++) {
 509       ret=(fabs(coord1[i]-coord2[i])<1e-15);
 510     }
 511   }
 512   if(ret) {
 513     const CONNECTIVITY* CNV = other.getConnectivityptr();
 514     if( (!CNV && _connectivity) || (CNV && !(_connectivity)) ) {
 515       return false;
 516     }
 517     if(_connectivity) {
 518       return _connectivity->deepCompare(*other._connectivity);
 519     }
 520   }
 521   return ret;
 522 }
 523 /*!
 524 @}
 525 */
 526
 527 /*!
 528  * \brief print my contents
 529  */
 530 ostream & ::MEDMEM::operator<<(ostream &os, const MESH &myMesh)
 531 {
 532   myMesh.printMySelf(os);
 533   return os;
 534 }
 535 void MESH::printMySelf(ostream &os) const
 536 {
 537   const MESH &myMesh = *this;
 538   int spacedimension = myMesh.getSpaceDimension();
 539   int meshdimension  = myMesh.getMeshDimension();
 540   int numberofnodes  = myMesh.getNumberOfNodes();
 541
 542   if ( spacedimension == MED_INVALID ) {
 543     os << " Empty mesh ...";
 544     return;
 545   }
 546
 547   os << "Space Dimension : " << spacedimension << endl << endl;
 548
 549   os << "Mesh Dimension : " << meshdimension << endl << endl;
 550
 551   if(myMesh.getCoordinateptr()) {
 552     const double * coordinates = myMesh.getCoordinates(MED_FULL_INTERLACE);
 553
 554     os << "SHOW NODES COORDINATES : " << endl;
 555     os << "Name :" << endl;
 556     const string * coordinatesnames = myMesh.getCoordinatesNames();
 557     if ( coordinatesnames ) {
 558       for(int i=0; i<spacedimension ; i++)
 559       {
 560         os << " - " << coordinatesnames[i] << endl;
 561       }
 562     }
 563     os << "Unit :" << endl;
 564     const string * coordinatesunits = myMesh.getCoordinatesUnits();
 565     if ( coordinatesunits ) {
 566       for(int i=0; i<spacedimension ; i++)
 567       {
 568         os << " - " << coordinatesunits[i] << endl;
 569       }
 570     }
 571     for(int i=0; i<numberofnodes ; i++)
 572     {
 573       os << "Nodes " << i+1 << " : ";
 574       for (int j=0; j<spacedimension ; j++)
 575         os << coordinates[i*spacedimension+j] << " ";
 576       os << endl;
 577     }
 578   }
 579
 580   if(myMesh.getConnectivityptr()) {
 581     os << endl << "SHOW CONNECTIVITY  :" << endl;
 582     os << *myMesh._connectivity << endl;
 583   }
 584
 585   medEntityMesh entity;
 586   os << endl << "SHOW FAMILIES :" << endl << endl;
 587   for (int k=1; k<=4; k++)
 588     {
 589       if (k==1) entity = MED_NODE;
 590       if (k==2) entity = MED_CELL;
 591       if (k==3) entity = MED_FACE;
 592       if (k==4) entity = MED_EDGE;
 593       int numberoffamilies = myMesh.getNumberOfFamilies(entity);
 594       os << "NumberOfFamilies on "<< entNames[entity] <<" : "<<numberoffamilies<<endl;
 595       using namespace MED_EN;
 596       for (int i=1; i<numberoffamilies+1;i++)
 597         {
 598           os << * myMesh.getFamily(entity,i) << endl;
 599         }
 600     }
 601
 602   os << endl << "SHOW GROUPS :" << endl << endl;
 603   for (int k=1; k<=4; k++)
 604     {
 605       if (k==1) entity = MED_NODE;
 606       if (k==2) entity = MED_CELL;
 607       if (k==3) entity = MED_FACE;
 608       if (k==4) entity = MED_EDGE;
 609       int numberofgroups = myMesh.getNumberOfGroups(entity);
 610       os << "NumberOfGroups on "<< entNames[entity] <<" : "<<numberofgroups<<endl;
 611       using namespace MED_EN;
 612       for (int i=1; i<numberofgroups+1;i++)
 613         {
 614           os << * myMesh.getGroup(entity,i) << endl;
 615         }
 616     }
 617 }
 618
 619 /*!
 620   Get global number of element which have same connectivity than connectivity argument.
 621
 622   It do not take care of connectivity order (3,4,7,10 is same as 7,3,10,4).
 623
 624   Return -1 if not found.
 625 */
 626 int MESH::getElementNumber(MED_EN::medConnectivity ConnectivityType,
 627                                                                                                          MED_EN::medEntityMesh Entity,
 628                                                                                                          MED_EN::medGeometryElement Type,
 629                                                                                                          int * connectivity) const
 630 {
 631   const char* LOC="MESH::getElementNumber " ;
 632   BEGIN_OF_MED(LOC) ;
 633
 634   int numberOfValue ; // size of connectivity array
 635   CELLMODEL myType(Type) ;
 636   if (ConnectivityType==MED_DESCENDING)
 637     numberOfValue = myType.getNumberOfConstituents(1) ; // faces (3D) or edges (2D)
 638   else
 639     numberOfValue = myType.getNumberOfNodes() ; // nodes
 640
 641   const int * myReverseConnectivityValue = getReverseConnectivity(ConnectivityType,Entity) ;
 642   const int * myReverseConnectivityIndex = getReverseConnectivityIndex(ConnectivityType,Entity) ;
 643
 644   // First node or face/edge
 645   int indexBegin = myReverseConnectivityIndex[connectivity[0]-1] ;
 646   int indexEnd = myReverseConnectivityIndex[connectivity[0]] ;
 647
 648   // list to put cells numbers
 649   list<int> cellsList ;
 650   list<int>::iterator itList ;
 651   for (int i=indexBegin; i<indexEnd; i++)
 652     cellsList.push_back(myReverseConnectivityValue[i-1]) ;
 653
 654   // Foreach node or face/edge in cell, we search cells which are in common.
 655   // At the end, we must have only one !
 656
 657   for (int i=1; i<numberOfValue; i++) {
 658     int connectivity_i = connectivity[i] ;
 659     for (itList=cellsList.begin();itList!=cellsList.end();itList++) {
 660       bool find = false ;
 661       for (int j=myReverseConnectivityIndex[connectivity_i-1]; j<myReverseConnectivityIndex[connectivity_i]; j++) {
 662         if ((*itList)==myReverseConnectivityValue[j-1]) {
 663           find=true ;
 664           break ;
 665         }
 666       }
 667       if (!find) {
 668         itList=cellsList.erase(itList);
 669         itList--; // well : rigth if itList = cellsList.begin() ??
 670       }
 671     }
 672   }
 673
 674   if (cellsList.size()>1) // we have more than one cell
 675     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Found more than one element !")) ;
 676
 677   if (cellsList.size()==0)
 678     return -1;
 679
 680   END_OF_MED(LOC);
 681
 682   return cellsList.front() ;
 683 }
 684
 685 /*!
 686 \addtogroup MESH_advanced
 687 @{
 688 The methods described in this section are algorithms that perform a computation
 689 and return a result in the form of a SUPPORT or a FIELD to the user. For large meshes,
 690 as the returned information is not stored in the mesh but is rather computed, the
 691 computation time can be large.
 692 */
 693
 694 /*!
 695   Returns a support which reference all elements on the boundary of mesh.
 696   For a d-dimensional mesh, a boundary element is defined as a d-1 dimension
 697   element that is referenced by only one element in the full descending connectivity.
 698
 699   This method can also return the list of nodes that belong to the boundary elements.
 700
 701   WARNING: This method can recalculate descending connectivity from partial to full form,
 702   so that partial SUPPORT on d-1 dimension elements becomes invalid.
 703
 704   \param Entity entity on which the boundary is desired. It has to be either \a MED_NODE or the
 705   d-1 dimension entity type (MED_FACE in 3D, MED_EDGE in 2D).
 706 */
 707 SUPPORT * MESH::getBoundaryElements(MED_EN::medEntityMesh Entity)
 708   throw (MEDEXCEPTION)
 709 {
 710   const char * LOC = "MESH::getBoundaryElements : " ;
 711   BEGIN_OF_MED(LOC) ;
 712   // some test :
 713   // actually we could only get face (in 3D) and edge (in 2D)
 714   medEntityMesh entityToParse=Entity;
 715   if(_spaceDimension == 3)
 716     if (Entity != MED_FACE)
 717       if(Entity==MED_NODE)
 718         entityToParse=MED_FACE;
 719       else
 720         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<Entity<<" !"));
 721   if(_spaceDimension == 2)
 722     if(Entity != MED_EDGE)
 723       if(Entity==MED_NODE)
 724         entityToParse=MED_EDGE;
 725       else
 726         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<Entity<<" !"));
 727
 728   // assure that descending connectivity is full
 729   if ( !_connectivity )
 730     throw MEDEXCEPTION("MESH::getgetBoundaryElements() : no connectivity defined in MESH !");
 731   _connectivity->calculateFullDescendingConnectivity(MED_CELL);
 732
 733   const int * myConnectivityValue = getReverseConnectivity(MED_DESCENDING) ;
 734   const int * myConnectivityIndex = getReverseConnectivityIndex(MED_DESCENDING) ;
 735   int numberOf = getNumberOfElementsWithPoly(entityToParse,MED_ALL_ELEMENTS) ;
 736   list<int> myElementsList;
 737
 738   for (int i=0 ; i<numberOf; i++)
 739     if (myConnectivityValue[myConnectivityIndex[i]] == 0) {
 740       myElementsList.push_back(i+1);
 741     }
 742   if ( myElementsList.empty() && numberOf != 0 )
 743     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No boundary elements found by reverse descending connectivity for entity "<<Entity<<" !"));
 744
 745   if(Entity==MED_NODE)
 746     return buildSupportOnNodeFromElementList(myElementsList,entityToParse);
 747   else
 748     return buildSupportOnElementsFromElementList(myElementsList,entityToParse);
 749 }
 750 /*!
 751 @}
 752 */
 753
 754 /*!
 755   Method return a reference on a support define on all the element of an entity.
 756 */
 757
 758 SUPPORT * MESH::getSupportOnAll(medEntityMesh entity)
 759   throw(MEDEXCEPTION)
 760 {
 761   const char * LOC = "MESH::getSupportOnAll : " ;
 762   BEGIN_OF_MED(LOC) ;
 763   if(entity == MED_ALL_ENTITIES)
 764     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined on entity MED_ALL_ENTITIES !"));
 765
 766   map<medEntityMesh,SUPPORT*>::const_iterator it =  _entitySupport.find(entity);
 767
 768   // find support and return is if exists
 769   if(it != _entitySupport.end())
 770     return (*it).second;
 771   else{
 772
 773     //build, store and return support
 774     string aSuppName = "SupportOnAll_"+entNames[entity];
 775     SUPPORT * aSupport = new SUPPORT((MESH *)this,aSuppName,entity);
 776
 777     _entitySupport.insert(make_pair(entity,aSupport));
 778     return aSupport;
 779   }
 780 }
 781
 782
 783 /*!
 784   Method that do the same thing as buildSupportOnNodeFromElementList except that a SUPPORT is not created.
 785  */
 786 void MESH::fillSupportOnNodeFromElementList(const list<int>& listOfElt, SUPPORT *supportToFill) const throw (MEDEXCEPTION)
 787 {
 788   MED_EN::medEntityMesh entity=supportToFill->getEntity();
 789   supportToFill->setAll(false);
 790   supportToFill->setMeshDirectly((MESH *)this);
 791
 792   int i;
 793   set<int> nodes;
 794   if ( entity == MED_NODE ) {
 795     supportToFill->fillFromNodeList(listOfElt);
 796   }
 797   else {
 798     for(list<int>::const_iterator iter=listOfElt.begin();iter!=listOfElt.end();iter++)
 799     {
 800       int lgth;
 801       const int *conn=_connectivity->getConnectivityOfAnElementWithPoly(MED_NODAL,entity,*iter,lgth);
 802       for(i=0;i<lgth;i++)
 803         nodes.insert(conn[i]);
 804     }
 805     list<int> nodesList;
 806     for(set<int>::iterator iter2=nodes.begin();iter2!=nodes.end();iter2++)
 807       nodesList.push_back(*iter2);
 808     supportToFill->fillFromNodeList(nodesList);
 809   }
 810 }
 811
 812 /*!
 813   Method created to factorize code. This method creates a new support on NODE (to deallocate) containing all the nodes id contained in elements 'listOfElt' of
 814   entity 'entity'.
 815  */
 816 SUPPORT *MESH::buildSupportOnNodeFromElementList(const list<int>& listOfElt,MED_EN::medEntityMesh entity) const throw (MEDEXCEPTION)
 817 {
 818   SUPPORT * mySupport = new SUPPORT((MESH *)this,"Boundary",entity);
 819   fillSupportOnNodeFromElementList(listOfElt,mySupport);
 820   return mySupport;
 821 }
 822
 823 /*!
 824   Method created to factorize code. This method creates a new support on entity 'entity' (to deallocate) containing all the entities contained in
 825   elements 'listOfElt' of entity 'entity'.
 826  */
 827 SUPPORT *MESH::buildSupportOnElementsFromElementList(const list<int>& listOfElt, MED_EN::medEntityMesh entity) const throw (MEDEXCEPTION)
 828 {
 829   const char* LOC = "MESH::buildSupportOnElementsFromElementList : ";
 830   BEGIN_OF_MED(LOC);
 831   SUPPORT *mySupport=new SUPPORT((MESH *)this,"Boundary",entity);
 832   mySupport->fillFromElementList(listOfElt);
 833   END_OF_MED(LOC);
 834   return mySupport ;
 835 }
 836
 837 /*!
 838 \addtogroup MESH_advanced
 839 @{
 840 */
 841 /*! Retrieves the volume of all the elements contained in \a Support. This method returns
 842 a FIELD structure based on this support. It only works on MED_CELL for 3D meshes.
 843 */
 844 FIELD<double, FullInterlace>* MESH::getVolume(const SUPPORT *Support) const throw (MEDEXCEPTION)
 845 {
 846   const char * LOC = "MESH::getVolume(SUPPORT*) : ";
 847   BEGIN_OF_MED(LOC);
 848
 849   // Support must be on 3D elements
 850
 851   // Make sure that the MESH class is the same as the MESH class attribut
 852   // in the class Support
 853   MESH* myMesh = Support->getMesh();
 854   if (this != myMesh)
 855     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
 856
 857   int dim_space = getSpaceDimension();
 858   medEntityMesh support_entity = Support->getEntity();
 859   bool onAll = Support->isOnAllElements();
 860
 861   int nb_type, length_values;
 862   const medGeometryElement* types;
 863   int nb_entity_type;
 864   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
 865   const int* global_connectivity;
 866   nb_type = Support->getNumberOfTypes();
 867   length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
 868   types = Support->getTypes();
 869   int index;
 870   FIELD<double, FullInterlace>* Volume =
 871     new FIELD<double, FullInterlace>(Support,1);
 872   //  double *volume = new double [length_values];
 873   Volume->setName("VOLUME");
 874   Volume->setDescription("cells volume");
 875   Volume->setComponentName(1,"volume");
 876   Volume->setComponentDescription(1,"desc-comp");
 877
 878   /*  string MEDComponentUnit(MED_TAILLE_PNOM,' ');*/
 879
 880   string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1]+"*"+myMesh->getCoordinatesUnits()[2];
 881
 882   Volume->setMEDComponentUnit(1,MEDComponentUnit);
 883
 884   Volume->setIterationNumber(0);
 885   Volume->setOrderNumber(0);
 886   Volume->setTime(0.0);
 887
 888   //const double *volume = Volume->getValue(MED_FULL_INTERLACE);
 889   typedef  MEDMEM_ArrayInterface<double,FullInterlace,NoGauss>::Array ArrayNoGauss;
 890   ArrayNoGauss  *volume = Volume->getArrayNoGauss();
 891
 892   index = 1;
 893   const double * coord = getCoordinates(MED_FULL_INTERLACE);
 894
 895   for (int i=0;i<nb_type;i++)
 896     {
 897       medGeometryElement type = types[i] ;
 898       double xvolume;
 899       nb_entity_type = Support->getNumberOfElements(type);
 900       if(type != MED_EN::MED_POLYHEDRA)
 901         {
 902           if (onAll)
 903             {
 904               global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
 905             }
 906           else
 907             {
 908               const int * supp_number = Support->getNumber(type);
 909               const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
 910               const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
 911               int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
 912
 913               for (int k_type = 0; k_type<nb_entity_type; k_type++) {
 914                 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
 915                   global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
 916                 }
 917               }
 918               global_connectivity = global_connectivity_tmp ;
 919             }
 920         }
 921
 922       switch (type)
 923         {
 924         case MED_TETRA4 : case MED_TETRA10 :
 925           {
 926             for (int tetra=0;tetra<nb_entity_type;tetra++)
 927               {
 928                 int tetra_index = (type%100)*tetra;
 929                 int N1 = global_connectivity[tetra_index]-1;
 930                 int N2 = global_connectivity[tetra_index+1]-1;
 931                 int N3 = global_connectivity[tetra_index+2]-1;
 932                 int N4 = global_connectivity[tetra_index+3]-1;
 933                 xvolume=INTERP_KERNEL::calculateVolumeForTetra(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4);
 934                 volume->setIJ(index,1,xvolume) ;
 935                 index++;
 936               }
 937             break;
 938           }
 939         case MED_PYRA5 : case MED_PYRA13 :
 940           {
 941             for (int pyra=0;pyra<nb_entity_type;pyra++)
 942               {
 943                 int pyra_index = (type%100)*pyra;
 944                 int N1 = global_connectivity[pyra_index]-1;
 945                 int N2 = global_connectivity[pyra_index+1]-1;
 946                 int N3 = global_connectivity[pyra_index+2]-1;
 947                 int N4 = global_connectivity[pyra_index+3]-1;
 948                 int N5 = global_connectivity[pyra_index+4]-1;
 949                 xvolume=INTERP_KERNEL::calculateVolumeForPyra(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4,coord+dim_space*N5);
 950                 volume->setIJ(index,1,xvolume) ;
 951                 index++;
 952               }
 953             break;
 954           }
 955         case MED_PENTA6 : case MED_PENTA15 :
 956           {
 957             for (int penta=0;penta<nb_entity_type;penta++)
 958               {
 959                 int penta_index = (type%100)*penta;
 960                 int N1 = global_connectivity[penta_index]-1;
 961                 int N2 = global_connectivity[penta_index+1]-1;
 962                 int N3 = global_connectivity[penta_index+2]-1;
 963                 int N4 = global_connectivity[penta_index+3]-1;
 964                 int N5 = global_connectivity[penta_index+4]-1;
 965                 int N6 = global_connectivity[penta_index+5]-1;
 966                 xvolume=INTERP_KERNEL::calculateVolumeForPenta(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4,coord+dim_space*N5,coord+dim_space*N6);
 967                 volume->setIJ(index,1,xvolume) ;
 968                 index++;
 969               }
 970             break;
 971           }
 972         case MED_HEXA8 : case MED_HEXA20 :
 973           {
 974             for (int hexa=0;hexa<nb_entity_type;hexa++)
 975               {
 976                 int hexa_index = (type%100)*hexa;
 977
 978                 int N1 = global_connectivity[hexa_index]-1;
 979                 int N2 = global_connectivity[hexa_index+1]-1;
 980                 int N3 = global_connectivity[hexa_index+2]-1;
 981                 int N4 = global_connectivity[hexa_index+3]-1;
 982                 int N5 = global_connectivity[hexa_index+4]-1;
 983                 int N6 = global_connectivity[hexa_index+5]-1;
 984                 int N7 = global_connectivity[hexa_index+6]-1;
 985                 int N8 = global_connectivity[hexa_index+7]-1;
 986                 xvolume=INTERP_KERNEL::calculateVolumeForHexa(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4,coord+dim_space*N5,coord+dim_space*N6,coord+dim_space*N7,coord+dim_space*N8);
 987                 volume->setIJ(index,1,xvolume) ;
 988                 index++;
 989               }
 990             break;
 991           }
 992         case MED_POLYHEDRA:
 993           {
 994             double bary[3];
 995             if(onAll)
 996               {
 997                 for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
 998                   {
 999                     int lgthNodes,iPts,iFaces,iPtsInFace;
1000                     int offsetWithClassicType=getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1001                     int *nodes=_connectivity->getNodesOfPolyhedron(offsetWithClassicType+polyhs+1,lgthNodes);
1002                     int nbOfFaces,*nbOfNodesPerFaces;
1003                     int **nodes1=_connectivity->getNodesPerFaceOfPolyhedron(offsetWithClassicType+polyhs+1,nbOfFaces,nbOfNodesPerFaces);
1004                     double **pts=new double * [lgthNodes];
1005                     double ***pts1=new double ** [nbOfFaces];
1006                     for(iPts=0;iPts<lgthNodes;iPts++)
1007                       pts[iPts]=(double *)(coord+3*(nodes[iPts]-1));
1008                     for(iFaces=0;iFaces<nbOfFaces;iFaces++)
1009                       {
1010                         pts1[iFaces]=new double* [nbOfNodesPerFaces[iFaces]];
1011                         for(iPtsInFace=0;iPtsInFace<nbOfNodesPerFaces[iFaces];iPtsInFace++)
1012                           pts1[iFaces][iPtsInFace]=(double *)(coord+3*(nodes1[iFaces][iPtsInFace]-1));
1013                       }
1014                     delete [] nodes1;
1015                     INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,lgthNodes,3,bary);
1016                     delete [] nodes;
1017                     delete [] pts;
1018                     xvolume=INTERP_KERNEL::calculateVolumeForPolyh((const double ***)pts1,nbOfNodesPerFaces,nbOfFaces,bary);
1019                     delete [] nbOfNodesPerFaces;
1020                     for(iFaces=0;iFaces<nbOfFaces;iFaces++)
1021                         delete [] pts1[iFaces];
1022                     delete [] pts1;
1023                     volume->setIJ(index,1,xvolume) ;
1024                     index++;
1025                   }
1026               }
1027             else
1028               {
1029                 const int * supp_number = Support->getNumber(MED_EN::MED_POLYHEDRA);
1030                 for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
1031                   {
1032                     int lgthNodes,iPts,iFaces,iPtsInFace;
1033                     int *nodes=_connectivity->getNodesOfPolyhedron(supp_number[polyhs],lgthNodes);
1034                     int nbOfFaces,*nbOfNodesPerFaces;
1035                     int **nodes1=_connectivity->getNodesPerFaceOfPolyhedron(supp_number[polyhs],nbOfFaces,nbOfNodesPerFaces);
1036                     double **pts=new double * [lgthNodes];
1037                     double ***pts1=new double ** [nbOfFaces];
1038                     for(iPts=0;iPts<lgthNodes;iPts++)
1039                       pts[iPts]=(double *)(coord+3*(nodes[iPts]-1));
1040                     for(iFaces=0;iFaces<nbOfFaces;iFaces++)
1041                       {
1042                         pts1[iFaces]=new double* [nbOfNodesPerFaces[iFaces]];
1043                         for(iPtsInFace=0;iPtsInFace<nbOfNodesPerFaces[iFaces];iPtsInFace++)
1044                           pts1[iFaces][iPtsInFace]=(double *)(coord+3*(nodes1[iFaces][iPtsInFace]-1));
1045                       }
1046                     delete [] nodes1;
1047                     INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,lgthNodes,3,bary);
1048                     delete [] nodes;
1049                     delete [] pts;
1050                     xvolume=INTERP_KERNEL::calculateVolumeForPolyh((const double ***)pts1,nbOfNodesPerFaces,nbOfFaces,bary);
1051                     delete [] nbOfNodesPerFaces;
1052                     for(iFaces=0;iFaces<nbOfFaces;iFaces++)
1053                         delete [] pts1[iFaces];
1054                     delete [] pts1;
1055                     volume->setIJ(index,1,xvolume) ;
1056                     index++;
1057                   }
1058               }
1059             break;
1060           }
1061         default :
1062           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Volumes on it !"));
1063           break;
1064         }
1065
1066       if (!onAll && type!=MED_EN::MED_POLYHEDRA)
1067         delete [] global_connectivity ;
1068     }
1069
1070   return Volume;
1071 }
1072 /*! Retrieves the area of all the elements contained in \a Support. This method returns
1073 a FIELD structure based on this support. It only works on MED_CELL for 2D meshes or MED_FACE
1074 for 3D meshes.
1075 */
1076
1077 FIELD<double, FullInterlace>* MESH::getArea(const SUPPORT * Support) const throw (MEDEXCEPTION)
1078 {
1079   const char * LOC = "MESH::getArea(SUPPORT*) : ";
1080   BEGIN_OF_MED(LOC);
1081
1082   // Support must be on 2D elements
1083
1084   // Make sure that the MESH class is the same as the MESH class attribut
1085   // in the class Support
1086   MESH* myMesh = Support->getMesh();
1087   if (this != myMesh)
1088     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1089
1090   int dim_space = getSpaceDimension();
1091   medEntityMesh support_entity = Support->getEntity();
1092   bool onAll = Support->isOnAllElements();
1093
1094   int nb_type, length_values;
1095   const medGeometryElement* types;
1096   int nb_entity_type;
1097   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1098   const int* global_connectivity;
1099
1100   nb_type = Support->getNumberOfTypes();
1101   length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1102   types = Support->getTypes();
1103
1104   int index;
1105   FIELD<double, FullInterlace>* Area;
1106
1107   Area = new FIELD<double, FullInterlace>(Support,1);
1108   Area->setName("AREA");
1109   Area->setDescription("cells or faces area");
1110
1111   Area->setComponentName(1,"area");
1112   Area->setComponentDescription(1,"desc-comp");
1113
1114   /*  string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1];*/
1115
1116   string MEDComponentUnit="trucmuch";
1117
1118   Area->setMEDComponentUnit(1,MEDComponentUnit);
1119
1120   Area->setIterationNumber(0);
1121   Area->setOrderNumber(0);
1122   Area->setTime(0.0);
1123
1124   double *area = (double *)Area->getValue();
1125
1126   const double * coord = getCoordinates(MED_FULL_INTERLACE);
1127   index = 0;
1128
1129   for (int i=0;i<nb_type;i++)
1130     {
1131       medGeometryElement type = types[i] ;
1132       nb_entity_type = Support->getNumberOfElements(type);
1133       const int *global_connectivityIndex = 0;
1134       if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
1135         {
1136           global_connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1137           if (onAll)
1138             {
1139               global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1140             }
1141           else
1142             {
1143               const int * supp_number = Support->getNumber(type);
1144               const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1145               int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1146
1147               for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1148                 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1149                   global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[global_connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1150                 }
1151               }
1152               global_connectivity = global_connectivity_tmp ;
1153             }
1154         }
1155       switch (type)
1156         {
1157         case MED_TRIA3 : case MED_TRIA6 :
1158           {
1159             for (int tria=0;tria<nb_entity_type;tria++)
1160               {
1161                 int tria_index = (type%100)*tria;
1162
1163                 int N1 = global_connectivity[tria_index]-1;
1164                 int N2 = global_connectivity[tria_index+1]-1;
1165                 int N3 = global_connectivity[tria_index+2]-1;
1166
1167                 area[index]=INTERP_KERNEL::calculateAreaForTria(coord+(dim_space*N1),
1168                                                    coord+(dim_space*N2),
1169                                                    coord+(dim_space*N3),dim_space);
1170                 index++;
1171               }
1172             break;
1173           }
1174         case MED_QUAD4 : case MED_QUAD8 :
1175           {
1176             for (int quad=0;quad<nb_entity_type;quad++)
1177               {
1178                 int quad_index = (type%100)*quad;
1179
1180                 int N1 = global_connectivity[quad_index]-1;
1181                 int N2 = global_connectivity[quad_index+1]-1;
1182                 int N3 = global_connectivity[quad_index+2]-1;
1183                 int N4 = global_connectivity[quad_index+3]-1;
1184
1185                 area[index]=INTERP_KERNEL::calculateAreaForQuad(coord+dim_space*N1,
1186                                                    coord+dim_space*N2,
1187                                                    coord+dim_space*N3,
1188                                                    coord+dim_space*N4,dim_space);
1189                 index++;
1190               }
1191             break;
1192           }
1193         case MED_POLYGON :
1194           {
1195             if(onAll)
1196               {
1197                 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
1198                 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
1199                 for (int polygs=0;polygs<nb_entity_type;polygs++)
1200                   {
1201                     int size=connectivity_index[polygs+1]-connectivity_index[polygs];
1202                     double **pts=new double * [size];
1203                     for(int iPts=0;iPts<size;iPts++)
1204                       pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[polygs]+iPts-1]-1));
1205                     area[index] = INTERP_KERNEL::calculateAreaForPolyg((const double **)pts,size,dim_space);
1206                     delete [] pts;
1207                     index++;
1208                   }
1209               }
1210             else
1211               {
1212                 const int * supp_number = Support->getNumber(MED_EN::MED_POLYGON);
1213                 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
1214                 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
1215                 int offsetWithClassicType=getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1216                 for (int polygs=0;polygs<nb_entity_type;polygs++)
1217                   {
1218                     int size=connectivity_index[supp_number[polygs]-offsetWithClassicType]-connectivity_index[supp_number[polygs]-offsetWithClassicType-1];
1219                     double **pts=new double * [size];
1220                     for(int iPts=0;iPts<size;iPts++)
1221                       pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[supp_number[polygs]-offsetWithClassicType-1]+iPts-1]-1));
1222                     area[index]=INTERP_KERNEL::calculateAreaForPolyg((const double **)pts,size,dim_space);
1223                     delete [] pts;
1224                     index++;
1225                   }
1226               }
1227             break;
1228           }
1229         default :
1230           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Areas on it !"));
1231           break;
1232         }
1233
1234       if (!onAll)
1235         if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
1236           delete [] global_connectivity ;
1237     }
1238   return Area;
1239 }
1240 /*! Retrieves the length of all the elements contained in \a Support. This method returns
1241 a FIELD structure based on this support. It only works on MED_EDGE supports.
1242 */
1243 FIELD<double, FullInterlace>* MESH::getLength(const SUPPORT * Support) const throw (MEDEXCEPTION)
1244 {
1245   const char * LOC = "MESH::getLength(SUPPORT*) : ";
1246   BEGIN_OF_MED(LOC);
1247
1248   // Support must be on 1D elements
1249
1250   // Make sure that the MESH class is the same as the MESH class attribut
1251   // in the class Support
1252   MESH* myMesh = Support->getMesh();
1253   if (this != myMesh)
1254     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1255
1256   int dim_space = getSpaceDimension();
1257   medEntityMesh support_entity = Support->getEntity();
1258   bool onAll = Support->isOnAllElements();
1259
1260   int nb_type, length_values;
1261   const medGeometryElement* types;
1262   int nb_entity_type;
1263   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1264   const int* global_connectivity;
1265
1266   nb_type = Support->getNumberOfTypes();
1267   length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1268   types = Support->getTypes();
1269
1270   int index;
1271   FIELD<double, FullInterlace>* Length;
1272
1273   Length = new FIELD<double, FullInterlace>(Support,1);
1274   //  double *length = new double [length_values];
1275
1276   //const double *length = Length->getValue(MED_FULL_INTERLACE);
1277   typedef  MEDMEM_ArrayInterface<double,FullInterlace,NoGauss>::Array ArrayNoGauss;
1278   ArrayNoGauss * length = Length->getArrayNoGauss();
1279
1280   const double * coord = getCoordinates(MED_FULL_INTERLACE);
1281   index = 1;
1282
1283   for (int i=0;i<nb_type;i++)
1284     {
1285       medGeometryElement type = types[i] ;
1286       double xlength;
1287
1288       if (onAll)
1289         {
1290           nb_entity_type = getNumberOfElements(support_entity,type);
1291           global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1292         }
1293       else
1294         {
1295           nb_entity_type = Support->getNumberOfElements(type);
1296           const int * supp_number = Support->getNumber(type);
1297           const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1298           const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1299           int* global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1300
1301           for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1302             for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1303               global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1304             }
1305           }
1306           global_connectivity = global_connectivity_tmp ;
1307         }
1308
1309       switch (type)
1310         {
1311         case MED_SEG2 : case MED_SEG3 :
1312           {
1313             for (int edge=0;edge<nb_entity_type;edge++)
1314               {
1315                 int edge_index = (type%100)*edge;
1316
1317                 int N1 = global_connectivity[edge_index]-1;
1318                 int N2 = global_connectivity[edge_index+1]-1;
1319
1320                 double x1 = coord[dim_space*N1];
1321                 double x2 = coord[dim_space*N2];
1322
1323                 double y1 = coord[dim_space*N1+1];
1324                 double y2 = coord[dim_space*N2+1];
1325
1326                 double z1, z2 ; z1 = 0.0 ; z2 = 0.0 ; if (dim_space==3) {
1327                   z1 = coord[dim_space*N1+2]; z2 = coord[dim_space*N2+2];}
1328
1329                 xlength =  sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) +
1330                                 (z1 - z2)*(z1 - z2));
1331
1332                 length->setIJ(index,1,xlength) ;
1333                 index++;
1334               }
1335             break;
1336           }
1337         default :
1338           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Lengths on it !"));
1339           break;
1340         }
1341
1342       if (!onAll) delete [] global_connectivity ;
1343     }
1344
1345   return Length;
1346 }
1347
1348 /*! Retrieves the normal for all elements contained in SUPPORT \a Support.
1349 The method is only functional for 2D supports for 3D meshes and 1D supports
1350 for 2D meshes. It returns
1351 a FIELD for which the number of components is equal to the dimension of the
1352 mesh and which represents coordinates of the vector normal to the element.
1353 The direction of the vector is undetermined.
1354 */
1355 FIELD<double, FullInterlace>* MESH::getNormal(const SUPPORT * Support) const throw (MEDEXCEPTION)
1356 {
1357   const char * LOC = "MESH::getNormal(SUPPORT*) : ";
1358   BEGIN_OF_MED(LOC);
1359
1360   // Support must be on 2D or 1D elements
1361
1362   // Make sure that the MESH class is the same as the MESH class attribut
1363   // in the class Support
1364   MESH* myMesh = Support->getMesh();
1365   if (this != myMesh)
1366     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : pointeur problem !"));
1367
1368   int dim_space = getSpaceDimension();
1369   int mesh_dim=getMeshDimension();
1370   medEntityMesh support_entity = Support->getEntity();
1371   bool onAll = Support->isOnAllElements();
1372
1373   if( support_entity!=MED_EDGE && (mesh_dim!=1 || support_entity!=MED_CELL) && ( mesh_dim!=2 || support_entity!=MED_CELL ) && ( mesh_dim!=3 || support_entity!=MED_FACE ))
1374     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"incompatible mesh dimension and entity"));
1375   int nb_type, length_values;
1376   const medGeometryElement* types;
1377   int nb_entity_type;
1378   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1379   const int* global_connectivity;
1380
1381   nb_type = Support->getNumberOfTypes();
1382   length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1383   types = Support->getTypes();
1384
1385   int index;
1386
1387   FIELD<double, FullInterlace>* Normal =
1388     new FIELD<double, FullInterlace>(Support,dim_space);
1389   Normal->setName("NORMAL");
1390   Normal->setDescription("cells or faces normal");
1391   for (int k=1;k<=dim_space;k++) {
1392     Normal->setComponentName(k,"normal");
1393     Normal->setComponentDescription(k,"desc-comp");
1394     Normal->setMEDComponentUnit(k,"unit");
1395   }
1396
1397   Normal->setIterationNumber(MED_NOPDT);
1398   Normal->setOrderNumber(MED_NONOR);
1399   Normal->setTime(0.0);
1400   double *normal = (double *)Normal->getValue();
1401
1402   const double * coord = getCoordinates(MED_FULL_INTERLACE);
1403   index = 0;
1404
1405   for (int i=0;i<nb_type;i++)
1406     {
1407       medGeometryElement type = types[i] ;
1408       nb_entity_type = Support->getNumberOfElements(type);
1409
1410       // Make sure we trying to get Normals on
1411       // TRIA3 or TRIA6 or QUAD4 QUAD8 (2D) cells on a 3D mesh
1412       // or on SEG2 or SEG3 (1D) cells on a 2D mesh
1413
1414       if ( (((type==MED_TRIA3) || (type==MED_TRIA6) ||
1415              (type==MED_QUAD4) || (type==MED_QUAD8) || (type==MED_POLYGON)) &&
1416             (dim_space != 3)) || (((type==MED_SEG2) || (type==MED_SEG3)) &&
1417                                   (dim_space != 2)) )
1418         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : dimension problem !"));
1419       if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
1420         {
1421           if (onAll)
1422             {
1423               global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1424             }
1425           else
1426             {
1427               const int * supp_number = Support->getNumber(type);
1428               const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1429               const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1430               int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1431
1432               for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1433                 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1434                   global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1435                 }
1436               }
1437
1438               global_connectivity = global_connectivity_tmp ;
1439             }
1440         }
1441
1442       switch (type)
1443         {
1444         case MED_TRIA3 : case MED_TRIA6 :
1445           {
1446             for (int tria=0;tria<nb_entity_type;tria++)
1447               {
1448                 int tria_index = (type%100)*tria;
1449                 int N1 = global_connectivity[tria_index]-1;
1450                 int N2 = global_connectivity[tria_index+1]-1;
1451                 int N3 = global_connectivity[tria_index+2]-1;
1452                 INTERP_KERNEL::calculateNormalForTria(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,normal+3*index);
1453                 index++;
1454               }
1455             break;
1456           }
1457         case MED_QUAD4 : case MED_QUAD8 :
1458           {
1459             for (int quad=0;quad<nb_entity_type;quad++)
1460               {
1461                 int quad_index = (type%100)*quad;
1462                 int N1 = global_connectivity[quad_index]-1;
1463                 int N2 = global_connectivity[quad_index+1]-1;
1464                 int N3 = global_connectivity[quad_index+2]-1;
1465                 int N4 = global_connectivity[quad_index+3]-1;
1466                 INTERP_KERNEL::calculateNormalForQuad(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4,normal+3*index);
1467                 index++;
1468               }
1469             break;
1470           }
1471         case MED_SEG2 : case MED_SEG3 :
1472           {
1473             double xnormal1, xnormal2;
1474             for (int edge=0;edge<nb_entity_type;edge++)
1475               {
1476                 int edge_index = (type%100)*edge;
1477
1478                 int N1 = global_connectivity[edge_index]-1;
1479                 int N2 = global_connectivity[edge_index+1]-1;
1480
1481                 double x1 = coord[dim_space*N1];
1482                 double x2 = coord[dim_space*N2];
1483
1484                 double y1 = coord[dim_space*N1+1];
1485                 double y2 = coord[dim_space*N2+1];
1486
1487                 xnormal1 = -(y2-y1);
1488                 xnormal2 = x2-x1;
1489
1490                 normal[2*index] = xnormal1 ;
1491                 normal[2*index+1] = xnormal2 ;
1492
1493                 index++;
1494               }
1495             break;
1496           }
1497         case MED_POLYGON :
1498           {
1499             if(onAll)
1500               {
1501                 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
1502                 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
1503                 for (int polygs=0;polygs<nb_entity_type;polygs++)
1504                   {
1505                     int size=connectivity_index[polygs+1]-connectivity_index[polygs];
1506                     double **pts=new double * [size];
1507                     for(int iPts=0;iPts<size;iPts++)
1508                       pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[polygs]+iPts-1])-1);
1509                     INTERP_KERNEL::calculateNormalForPolyg((const double **)pts,size,normal+3*index);
1510                     delete [] pts;
1511                     index++;
1512                   }
1513               }
1514             else
1515               {
1516                 const int * supp_number = Support->getNumber(MED_EN::MED_POLYGON);
1517                 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
1518                 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
1519                 int offsetWithClassicType=getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1520                 for (int polygs=0;polygs<nb_entity_type;polygs++)
1521                   {
1522                     int localPolygsNbP1=supp_number[polygs]-offsetWithClassicType;
1523                     int size=connectivity_index[localPolygsNbP1]-connectivity_index[localPolygsNbP1-1];
1524                     double **pts=new double * [size];
1525                     for(int iPts=0;iPts<size;iPts++)
1526                       pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[localPolygsNbP1-1]+iPts-1])-1);
1527                     INTERP_KERNEL::calculateNormalForPolyg((const double **)pts,size,normal+3*index);
1528                     delete [] pts;
1529                     index++;
1530                   }
1531               }
1532             break;
1533           }
1534         default :
1535           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Normals on it !"));
1536           break;
1537         }
1538       if (!onAll && type!=MED_EN::MED_POLYGON)
1539         delete [] global_connectivity ;
1540     }
1541   END_OF_MED(LOC);
1542
1543   return Normal;
1544 }
1545 /*!Returns the barycenter for each element in the support. The barycenter positions are returned
1546 as a field with a number of components equal to the mesh dimension.
1547 The barycenter computed by this method is actually the barycenter of the set of nodes of the elements, each having the same weight.
1548 */
1549 FIELD<double, FullInterlace>* MESH::getBarycenter(const SUPPORT * Support) const throw (MEDEXCEPTION)
1550 {
1551   const char * LOC = "MESH::getBarycenter(SUPPORT*) : ";
1552   MESH* myMesh = Support->getMesh();
1553   if (this != myMesh)
1554     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1555
1556   int dim_space = getSpaceDimension();
1557   medEntityMesh support_entity = Support->getEntity();
1558   bool onAll = Support->isOnAllElements();
1559
1560   int nb_type, length_values;
1561   const medGeometryElement* types;
1562   int nb_entity_type;
1563   const int* global_connectivity;
1564
1565   nb_type = Support->getNumberOfTypes();
1566   length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1567   types = Support->getTypes();
1568
1569   FIELD<double, FullInterlace>* Barycenter;
1570   Barycenter = new FIELD<double, FullInterlace>(Support,dim_space);
1571   Barycenter->setName("BARYCENTER");
1572   Barycenter->setDescription("cells or faces barycenter");
1573
1574   for (int k=0;k<dim_space;k++) {
1575     int kp1 = k + 1;
1576     Barycenter->setComponentName(kp1,myMesh->getCoordinatesNames()[k]);
1577     Barycenter->setComponentDescription(kp1,"desc-comp");
1578     Barycenter->setMEDComponentUnit(kp1,myMesh->getCoordinatesUnits()[k]);
1579   }
1580   Barycenter->setIterationNumber(0);
1581   Barycenter->setOrderNumber(0);
1582   Barycenter->setTime(0.0);
1583   double *barycenter=(double *)Barycenter->getValue();
1584   const double * coord = getCoordinates(MED_FULL_INTERLACE);
1585   int index=0;
1586   for (int i=0;i<nb_type;i++)
1587     {
1588       medGeometryElement type = types[i] ;
1589       nb_entity_type = Support->getNumberOfElements(type);
1590       if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA )
1591         {
1592           if (onAll)
1593             {
1594               global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1595             }
1596           else
1597             {
1598               const int * supp_number = Support->getNumber(type);
1599               const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1600               int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1601
1602               for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1603                 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1604                   const int *global_connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1605                   global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[global_connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1606                 }
1607               }
1608               global_connectivity = global_connectivity_tmp;
1609             }
1610         }
1611
1612       switch (type)
1613         {
1614         case MED_TETRA4 : case MED_TETRA10 :
1615           {
1616             for (int tetra =0;tetra<nb_entity_type;tetra++)
1617               {
1618                 int tetra_index = (type%100)*tetra;
1619
1620                 int N1 = global_connectivity[tetra_index]-1;
1621                 int N2 = global_connectivity[tetra_index+1]-1;
1622                 int N3 = global_connectivity[tetra_index+2]-1;
1623                 int N4 = global_connectivity[tetra_index+3]-1;
1624                 double *pts[4];
1625                 pts[0]=(double *)coord+dim_space*N1;
1626                 pts[1]=(double *)coord+dim_space*N2;
1627                 pts[2]=(double *)coord+dim_space*N3;
1628                 pts[3]=(double *)coord+dim_space*N4;
1629                 INTERP_KERNEL::calculateBarycenter<4,3>((const double **)pts,barycenter+3*index);
1630                 index++;
1631               }
1632             break;
1633           }
1634         case MED_PYRA5 : case MED_PYRA13 :
1635           {
1636             for (int pyra=0;pyra<nb_entity_type;pyra++)
1637               {
1638                 int pyra_index = (type%100)*pyra;
1639
1640                 int N1 = global_connectivity[pyra_index]-1;
1641                 int N2 = global_connectivity[pyra_index+1]-1;
1642                 int N3 = global_connectivity[pyra_index+2]-1;
1643                 int N4 = global_connectivity[pyra_index+3]-1;
1644                 int N5 = global_connectivity[pyra_index+4]-1;
1645                 double *pts[5];
1646                 pts[0]=(double *)coord+dim_space*N1;
1647                 pts[1]=(double *)coord+dim_space*N2;
1648                 pts[2]=(double *)coord+dim_space*N3;
1649                 pts[3]=(double *)coord+dim_space*N4;
1650                 pts[4]=(double *)coord+dim_space*N5;
1651                 INTERP_KERNEL::calculateBarycenter<5,3>((const double **)pts,barycenter+3*index);
1652                 index++;
1653               }
1654             break;
1655           }
1656         case MED_PENTA6 : case MED_PENTA15 :
1657           {
1658             for (int penta=0;penta<nb_entity_type;penta++)
1659               {
1660                 int penta_index = (type%100)*penta;
1661
1662                 int N1 = global_connectivity[penta_index]-1;
1663                 int N2 = global_connectivity[penta_index+1]-1;
1664                 int N3 = global_connectivity[penta_index+2]-1;
1665                 int N4 = global_connectivity[penta_index+3]-1;
1666                 int N5 = global_connectivity[penta_index+4]-1;
1667                 int N6 = global_connectivity[penta_index+5]-1;
1668                 double *pts[6];
1669                 pts[0]=(double *)coord+dim_space*N1;
1670                 pts[1]=(double *)coord+dim_space*N2;
1671                 pts[2]=(double *)coord+dim_space*N3;
1672                 pts[3]=(double *)coord+dim_space*N4;
1673                 pts[4]=(double *)coord+dim_space*N5;
1674                 pts[5]=(double *)coord+dim_space*N6;
1675                 INTERP_KERNEL::calculateBarycenter<6,3>((const double **)pts,barycenter+3*index);
1676                 index++;
1677               }
1678             break;
1679           }
1680         case MED_HEXA8 : case MED_HEXA20 :
1681           {
1682             for (int hexa=0;hexa<nb_entity_type;hexa++)
1683               {
1684                 int hexa_index = (type%100)*hexa;
1685
1686                 int N1 = global_connectivity[hexa_index]-1;
1687                 int N2 = global_connectivity[hexa_index+1]-1;
1688                 int N3 = global_connectivity[hexa_index+2]-1;
1689                 int N4 = global_connectivity[hexa_index+3]-1;
1690                 int N5 = global_connectivity[hexa_index+4]-1;
1691                 int N6 = global_connectivity[hexa_index+5]-1;
1692                 int N7 = global_connectivity[hexa_index+6]-1;
1693                 int N8 = global_connectivity[hexa_index+7]-1;
1694                 double *pts[8];
1695                 pts[0]=(double *)coord+dim_space*N1;
1696                 pts[1]=(double *)coord+dim_space*N2;
1697                 pts[2]=(double *)coord+dim_space*N3;
1698                 pts[3]=(double *)coord+dim_space*N4;
1699                 pts[4]=(double *)coord+dim_space*N5;
1700                 pts[5]=(double *)coord+dim_space*N6;
1701                 pts[6]=(double *)coord+dim_space*N7;
1702                 pts[7]=(double *)coord+dim_space*N8;
1703                 INTERP_KERNEL::calculateBarycenter<8,3>((const double **)pts,barycenter+3*index);
1704                 index++;
1705               }
1706             break;
1707           }
1708         case MED_TRIA3 : case MED_TRIA6 :
1709           {
1710             for (int tria=0;tria<nb_entity_type;tria++)
1711               {
1712                 int tria_index = (type%100)*tria;
1713                 int N1 = global_connectivity[tria_index]-1;
1714                 int N2 = global_connectivity[tria_index+1]-1;
1715                 int N3 = global_connectivity[tria_index+2]-1;
1716                 double *pts[3];
1717                 pts[0]=(double *)coord+dim_space*N1;
1718                 pts[1]=(double *)coord+dim_space*N2;
1719                 pts[2]=(double *)coord+dim_space*N3;
1720                 if (dim_space==2)
1721                   INTERP_KERNEL::calculateBarycenter<3,2>((const double **)pts,barycenter+2*index);
1722                 else
1723                   INTERP_KERNEL::calculateBarycenter<3,3>((const double **)pts,barycenter+3*index);
1724                 index++;
1725               }
1726             break;
1727           }
1728         case MED_QUAD4 : case MED_QUAD8 :
1729           {
1730             for (int quad=0;quad<nb_entity_type;quad++)
1731               {
1732                 int quad_index = (type%100)*quad;
1733                 int N1 = global_connectivity[quad_index]-1;
1734                 int N2 = global_connectivity[quad_index+1]-1;
1735                 int N3 = global_connectivity[quad_index+2]-1;
1736                 int N4 = global_connectivity[quad_index+3]-1;
1737                 double *pts[4];
1738                 pts[0]=(double *)coord+dim_space*N1;
1739                 pts[1]=(double *)coord+dim_space*N2;
1740                 pts[2]=(double *)coord+dim_space*N3;
1741                 pts[3]=(double *)coord+dim_space*N4;
1742                 if (dim_space==2)
1743                   INTERP_KERNEL::calculateBarycenter<4,2>((const double **)pts,barycenter+2*index);
1744                 else
1745                   INTERP_KERNEL::calculateBarycenter<4,3>((const double **)pts,barycenter+3*index);
1746                 index++;
1747               }
1748             break;
1749           }
1750         case MED_SEG2 : case MED_SEG3 :
1751           {
1752             for (int edge=0;edge<nb_entity_type;edge++)
1753               {
1754                 int edge_index = (type%100)*edge;
1755                 int N1 = global_connectivity[edge_index]-1;
1756                 int N2 = global_connectivity[edge_index+1]-1;
1757                 double *pts[2];
1758                 pts[0]=(double *)coord+dim_space*N1;
1759                 pts[1]=(double *)coord+dim_space*N2;
1760                 if (dim_space==2)
1761                   INTERP_KERNEL::calculateBarycenter<2,2>((const double **)pts,barycenter+2*index);
1762                 else
1763                   INTERP_KERNEL::calculateBarycenter<2,3>((const double **)pts,barycenter+3*index);
1764                 index++;
1765               }
1766             break;
1767           }
1768         case MED_POLYGON :
1769           {
1770             if(onAll)
1771               {
1772                 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
1773                 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
1774                 for (int polygs=0;polygs<nb_entity_type;polygs++)
1775                   {
1776                     int size=connectivity_index[polygs+1]-connectivity_index[polygs];
1777                     double **pts=new double * [size];
1778                     for(int iPts=0;iPts<size;iPts++)
1779                       pts[iPts]=(double *)coord+dim_space*(connectivity[connectivity_index[polygs]+iPts-1]-1);
1780                     INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,size,dim_space,barycenter+dim_space*index);
1781                     delete [] pts;
1782                   }
1783               }
1784             else
1785               {
1786                 const int * supp_number = Support->getNumber(MED_EN::MED_POLYGON);
1787                 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
1788                 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
1789                 int offsetWithClassicType=getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1790                 for (int polygs=0;polygs<nb_entity_type;polygs++)
1791                   {
1792                     int localPolygsNbP1=supp_number[polygs]-offsetWithClassicType;
1793                     int size=connectivity_index[localPolygsNbP1]-connectivity_index[localPolygsNbP1-1];
1794                     double **pts=new double * [size];
1795                     for(int iPts=0;iPts<size;iPts++)
1796                       pts[iPts]=(double *)coord+dim_space*(connectivity[connectivity_index[localPolygsNbP1-1]+iPts-1]-1);
1797                     INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,size,dim_space,barycenter+dim_space*index);
1798                     delete [] pts;
1799                   }
1800               }
1801             index++;
1802             break;
1803           }
1804         case MED_EN::MED_POLYHEDRA:
1805           {
1806             if(onAll)
1807               {
1808                 for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
1809                   {
1810                     int lgthNodes;
1811                     int offsetWithClassicType=getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1812                     int *nodes=_connectivity->getNodesOfPolyhedron(offsetWithClassicType+polyhs+1,lgthNodes);
1813                     double **pts=new double * [lgthNodes];
1814                     for(int iPts=0;iPts<lgthNodes;iPts++)
1815                       pts[iPts]=(double *)coord+3*(nodes[iPts]-1);
1816                     INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,lgthNodes,3,barycenter+3*index);
1817                     delete [] pts;
1818                     delete [] nodes;
1819                     index++;
1820                   }
1821               }
1822             else
1823               {
1824                 const int * supp_number = Support->getNumber(MED_EN::MED_POLYHEDRA);
1825                 for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
1826                   {
1827                     int lgthNodes;
1828                     int *nodes=_connectivity->getNodesOfPolyhedron(supp_number[polyhs],lgthNodes);
1829                     double **pts=new double * [lgthNodes];
1830                     for(int iPts=0;iPts<lgthNodes;iPts++)
1831                       pts[iPts]=(double *)coord+3*(nodes[iPts]-1);
1832                     INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,lgthNodes,3,barycenter+3*index);
1833                     delete [] pts;
1834                     delete [] nodes;
1835                     index++;
1836                   }
1837               }
1838             break;
1839           }
1840         default :
1841           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get a barycenter on it (in fact unknown type) !"));
1842           break;
1843         }
1844
1845       if (!onAll)
1846         if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
1847           delete [] global_connectivity;
1848     }
1849   //END_OF_MED();
1850   return Barycenter;
1851 }
1852 /*!
1853 @}
1854 */
1855
1856 bool MESH::isEmpty() const
1857 {
1858     bool notempty = _name != "NOT DEFINED"                || _coordinate != NULL           || _connectivity != NULL ||
1859                  _spaceDimension !=MED_INVALID || _meshDimension !=MED_INVALID  ||
1860                  _numberOfNodes !=MED_INVALID  || _groupNode.size() != 0   ||
1861                  _familyNode.size() != 0       || _groupCell.size() != 0   ||
1862                  _familyCell.size() != 0       || _groupFace.size() != 0   ||
1863                  _familyFace.size() != 0       || _groupEdge.size() != 0   ||
1864                  _familyEdge.size() != 0       || _isAGrid != 0 ;
1865     return !notempty;
1866 }
1867
1868 void MESH::read(int index)
1869 {
1870   const char * LOC = "MESH::read(int index=0) : ";
1871   BEGIN_OF_MED(LOC);
1872
1873   if (_drivers[index]) {
1874     _drivers[index]->open();
1875     _drivers[index]->read();
1876     _drivers[index]->close();
1877   }
1878   else
1879     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
1880                                      << "The index given is invalid, index must be between  0 and |"
1881                                      << _drivers.size()
1882                                      )
1883                           );
1884   END_OF_MED(LOC);
1885 }
1886
1887 /*!
1888 \addtogroup MESH_io
1889 @{
1890  */
1891 /*! Writes all the content of the MESH using driver referenced by the integer handle returned by a \a addDriver call.
1892
1893 Example :
1894 \verbatim
1895 //...
1896 // Attaching the driver to file "output.med", meshname "Mesh"
1897 int driver_handle = mesh.addDriver(MED_DRIVER, "output.med", "Mesh");
1898 // Writing the content of mesh to the file
1899 mesh.write(driver_handle);
1900 \endverbatim
1901 */
1902 void MESH::write(int index/*=0*/, const string & driverName/* = ""*/)
1903 {
1904   const char * LOC = "MESH::write(int index=0, const string & driverName = \"\") : ";
1905   BEGIN_OF_MED(LOC);
1906
1907   if ( _drivers[index] ) {
1908     _drivers[index]->open();
1909     if (driverName != "") _drivers[index]->setMeshName(driverName);
1910     _drivers[index]->write();
1911     _drivers[index]->close();
1912   }
1913   else
1914     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
1915                                      << "The index given is invalid, index must be between  0 and |"
1916                                      << _drivers.size()
1917                                      )
1918                           );
1919   END_OF_MED(LOC);
1920 }
1921 /*!
1922 @}
1923 */
1924
1925 /*!
1926 \addtogroup MESH_advanced
1927 @{
1928  */
1929
1930 /*!
1931 Retrieves the skin of support \a Support3D. This method is only available in 3D.
1932 On output, it returns a MED_FACE support with the skin of all elements contained in support.
1933 The skin is defined as the list of faces that are compnents of only one volume in the input
1934 support.
1935
1936 WARNING: This method can recalculate descending connectivity from partial to full form,
1937 so that partial SUPPORT on MED_FACE on this mesh becomes invalid.
1938  */
1939 SUPPORT * MESH::getSkin(const SUPPORT * Support3D) throw (MEDEXCEPTION)
1940 {
1941   const char * LOC = "MESH::getSkin : " ;
1942   BEGIN_OF_MED(LOC) ;
1943   // some test :
1944   if (this != Support3D->getMesh())
1945     throw MEDEXCEPTION(STRING(LOC) <<  "no compatibility between *this and SUPPORT::_mesh !");
1946   if (_meshDimension != 3 || Support3D->getEntity() != MED_CELL)
1947     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Defined on 3D cells only"));
1948
1949   // well, all rigth !
1950   SUPPORT * mySupport = new SUPPORT(this,"Skin",MED_FACE);
1951   mySupport->setAll(false);
1952
1953   list<int> myElementsList;
1954   int i,j, size = 0;
1955
1956   // assure that descending connectivity is full
1957   if ( !_connectivity )
1958     throw MEDEXCEPTION(STRING(LOC) << "no connectivity defined in MESH !");
1959   _connectivity->calculateFullDescendingConnectivity(MED_CELL);
1960   //calculateConnectivity(MED_FULL_INTERLACE, MED_DESCENDING, MED_CELL);
1961   if (Support3D->isOnAllElements())
1962   {
1963     const int* value = getReverseConnectivity(MED_DESCENDING);
1964     const int* index = getReverseConnectivityIndex(MED_DESCENDING);
1965
1966     int nbFaces = getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS);
1967     for (int i = 0; i < nbFaces; i++)
1968     {
1969       //a face is in skin if it is linked to one element or if one of the elements
1970       //it is linked to is "virtual"
1971       if ((index[i+1]-index[i]==1) || (value[index[i]-1]==0) || (value[index[i]]==0)) {
1972         myElementsList.push_back( i+1 );
1973         size++;
1974       }
1975     }
1976   }
1977   else
1978   {
1979     map<int,int> FaceNbEncounterNb;
1980
1981     int * myConnectivityValue = const_cast <int *>
1982       (getConnectivity(MED_FULL_INTERLACE, MED_DESCENDING,
1983                        MED_CELL, MED_ALL_ELEMENTS));
1984     int * myConnectivityIndex = const_cast <int *> (getConnectivityIndex(MED_DESCENDING, MED_CELL));
1985     int * myCellNbs = const_cast <int *> (Support3D->getnumber()->getValue());
1986     int nbCells = Support3D->getnumber()->getLength();
1987     for (i=0; i<nbCells; ++i)
1988     {
1989       int cellNb = myCellNbs [ i ];
1990       int faceFirst = myConnectivityIndex[ cellNb-1 ];
1991       int faceLast  = myConnectivityIndex[ cellNb ];
1992       for (j = faceFirst; j < faceLast; ++j)
1993       {
1994         int faceNb = abs( myConnectivityValue [ j-1 ] );
1995         //MESSAGE_MED( "Cell # " << i << " -- Face: " << faceNb);
1996         if (FaceNbEncounterNb.find( faceNb ) == FaceNbEncounterNb.end())
1997           FaceNbEncounterNb[ faceNb ] = 1;
1998         else
1999           FaceNbEncounterNb[ faceNb ] ++;
2000       }
2001     }
2002     map<int,int>::iterator FaceNbEncounterNbItr;
2003     for (FaceNbEncounterNbItr = FaceNbEncounterNb.begin();
2004          FaceNbEncounterNbItr != FaceNbEncounterNb.end();
2005          FaceNbEncounterNbItr ++)
2006       if ((*FaceNbEncounterNbItr).second == 1)
2007       {
2008         myElementsList.push_back( (*FaceNbEncounterNbItr).first) ;
2009         size++ ;
2010       }
2011   }
2012   // Well, we must know how many geometric type we have found
2013   int * myListArray = new int[size] ;
2014   int id = 0 ;
2015   list<int>::iterator myElementsListIt ;
2016   for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
2017     myListArray[id]=(*myElementsListIt) ;
2018     id ++ ;
2019   }
2020
2021   int numberOfGeometricType ;
2022   medGeometryElement* geometricType ;
2023   int * geometricTypeNumber ;
2024   int * numberOfEntities ;
2025   //  MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY() ;
2026   int * mySkyLineArrayIndex ;
2027
2028   int numberOfType = getNumberOfTypes(MED_FACE) ;
2029   if (numberOfType == 1) { // wonderfull : it's easy !
2030     numberOfGeometricType = 1 ;
2031     geometricType = new medGeometryElement[1] ;
2032     const medGeometryElement *  allType = getTypes(MED_FACE);
2033     geometricType[0] = allType[0] ;
2034     geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
2035     geometricTypeNumber[0] = 0 ;
2036     numberOfEntities = new int[1] ;
2037     numberOfEntities[0] = size ;
2038     mySkyLineArrayIndex = new int[2] ;
2039     mySkyLineArrayIndex[0]=1 ;
2040     mySkyLineArrayIndex[1]=1+size ;
2041   }
2042   else {// hemmm
2043     map<medGeometryElement,int> theType ;
2044     for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
2045       medGeometryElement myType = getElementType(MED_FACE,*myElementsListIt) ;
2046       if (theType.find(myType) != theType.end() )
2047         theType[myType]+=1 ;
2048       else
2049         theType[myType]=1 ;
2050     }
2051     numberOfGeometricType = theType.size() ;
2052     geometricType = new medGeometryElement[numberOfGeometricType] ;
2053     //const medGeometryElement *  allType = getTypes(MED_FACE); !! UNUSED VARIABLE !!
2054     geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
2055     numberOfEntities = new int[numberOfGeometricType] ;
2056     mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
2057     int index = 0 ;
2058     mySkyLineArrayIndex[0]=1 ;
2059     map<medGeometryElement,int>::iterator theTypeIt ;
2060     for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
2061       geometricType[index] = (*theTypeIt).first ;
2062       geometricTypeNumber[index] = 0 ;
2063       numberOfEntities[index] = (*theTypeIt).second ;
2064       mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfEntities[index] ;
2065       index++ ;
2066     }
2067   }
2068   //  mySkyLineArray->setMEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
2069   MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
2070
2071   mySupport->setNumberOfGeometricType(numberOfGeometricType) ;
2072   mySupport->setGeometricType(geometricType) ;
2073   //  mySupport->setGeometricTypeNumber(geometricTypeNumber) ;
2074   mySupport->setNumberOfElements(numberOfEntities) ;
2075   //mySupport->setTotalNumberOfElements(size) ;
2076   mySupport->setNumber(mySkyLineArray) ;
2077
2078   delete[] numberOfEntities;
2079   delete[] geometricTypeNumber;
2080   delete[] geometricType;
2081   delete[] mySkyLineArrayIndex;
2082   delete[] myListArray;
2083 //   delete mySkyLineArray;
2084
2085   END_OF_MED(LOC);
2086   return mySupport ;
2087
2088 }
2089
2090 /*!
2091   return a SUPPORT pointer on the union of all SUPPORTs in Supports.
2092   You should delete this pointer after use to avoid memory leaks.
2093 */
2094 SUPPORT * MESH::mergeSupports(const vector<SUPPORT *> Supports) throw (MEDEXCEPTION)
2095 {
2096   const char * LOC = "MESH:::mergeSupports(const vector<SUPPORT *> ) : " ;
2097   BEGIN_OF_MED(LOC) ;
2098
2099   SUPPORT * returnedSupport;
2100   string returnedSupportName;
2101   string returnedSupportDescription;
2102   char * returnedSupportNameChar;
2103   char * returnedSupportDescriptionChar;
2104   int size = Supports.size();
2105
2106   if (size == 0)
2107     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC) <<
2108                                      " mergeSupports() does't accept zero size vector"));
2109
2110   if (size == 1)
2111     {
2112       MESSAGE_MED(LOC <<" there is only one SUPPORT in the argument list, the method return a copy of this object !");
2113       SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2114
2115       returnedSupport = new SUPPORT(*obj);
2116
2117       int lenName = strlen((Supports[0]->getName()).c_str()) + 8 + 1;
2118       int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 8 + 1;
2119
2120       returnedSupportNameChar = new char[lenName];
2121       returnedSupportDescriptionChar = new char[lenDescription];
2122
2123       returnedSupportNameChar = strcpy(returnedSupportNameChar,"Copy of ");
2124       returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
2125       returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Copy of ");
2126       returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2127                                               (Supports[0]->getDescription()).c_str());
2128
2129       returnedSupportName = string(returnedSupportNameChar);
2130       returnedSupportDescription = string(returnedSupportDescriptionChar);
2131
2132       returnedSupport->setName(returnedSupportName);
2133       returnedSupport->setDescription(returnedSupportDescription);
2134
2135       delete [] returnedSupportNameChar;
2136       delete [] returnedSupportDescriptionChar;
2137     }
2138   else
2139     {
2140       SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2141       returnedSupport = new SUPPORT(*obj);
2142
2143       int lenName = strlen((Supports[0]->getName()).c_str()) + 9 + 1;
2144       int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 9 + 1;
2145
2146       for(int i = 1;i<size;i++)
2147         {
2148           obj = const_cast <SUPPORT *> (Supports[i]);
2149           returnedSupport->blending(obj);
2150
2151           if (i == (size-1))
2152             {
2153               lenName = lenName + 5 + strlen((Supports[i]->getName()).c_str());
2154               lenDescription = lenDescription + 5 +
2155                 strlen((Supports[i]->getDescription()).c_str());
2156             }
2157           else
2158             {
2159               lenName = lenName + 2 + strlen((Supports[i]->getName()).c_str());
2160               lenDescription = lenDescription + 2 +
2161                 strlen((Supports[i]->getDescription()).c_str());
2162             }
2163         }
2164
2165       returnedSupportNameChar = new char[lenName];
2166       returnedSupportDescriptionChar = new char[lenDescription];
2167
2168       returnedSupportNameChar = strcpy(returnedSupportNameChar,"Merge of ");
2169       returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Merge of ");
2170
2171       returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
2172       returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2173                                               (Supports[0]->getDescription()).c_str());
2174
2175       for(int i = 1;i<size;i++)
2176         {
2177           if (i == (size-1))
2178             {
2179               returnedSupportNameChar = strcat(returnedSupportNameChar," and ");
2180               returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar," and ");
2181
2182               returnedSupportNameChar = strcat(returnedSupportNameChar,
2183                                                (Supports[i]->getName()).c_str());
2184               returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2185                                                       (Supports[i]->getDescription()).c_str());
2186             }
2187           else
2188             {
2189               returnedSupportNameChar = strcat(returnedSupportNameChar,", ");
2190               returnedSupportNameChar = strcat(returnedSupportNameChar,
2191                                                (Supports[i]->getName()).c_str());
2192
2193               returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,", ");
2194               returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2195                                                       (Supports[i]->getDescription()).c_str());
2196             }
2197         }
2198
2199       returnedSupportName = string(returnedSupportNameChar);
2200       returnedSupport->setName(returnedSupportName);
2201
2202       returnedSupportDescription = string(returnedSupportDescriptionChar);
2203       returnedSupport->setDescription(returnedSupportDescription);
2204
2205       delete [] returnedSupportNameChar;
2206       delete [] returnedSupportDescriptionChar;
2207     }
2208
2209   END_OF_MED(LOC);
2210   return returnedSupport;
2211 }
2212
2213 /*!
2214   return a SUPPORT pointer on the intersection of all SUPPORTs in Supports.
2215   The (SUPPORT *) NULL pointer is returned if the intersection is empty.
2216   You should delete this pointer after use to avois memory leaks.
2217 */
2218 SUPPORT * MESH::intersectSupports(const vector<SUPPORT *> Supports) throw (MEDEXCEPTION)
2219 {
2220   const char* LOC = "MESH:::intersectSupports(const vector<SUPPORT *> ) : ";
2221   BEGIN_OF_MED(LOC);
2222
2223   SUPPORT * returnedSupport;
2224   string returnedSupportName;
2225   string returnedSupportDescription;
2226   char * returnedSupportNameChar;
2227   char * returnedSupportDescriptionChar;
2228   int size = Supports.size();
2229
2230   if (size == 1)
2231     {
2232       MESSAGE_MED(PREFIX_MED <<" there is only one SUPPORT in the argument list, the method return a copy of this object !");
2233       SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2234
2235       returnedSupport = new SUPPORT(*obj);
2236
2237       int lenName = strlen((Supports[0]->getName()).c_str()) + 8 + 1;
2238       int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 8 + 1;
2239
2240       returnedSupportNameChar = new char[lenName];
2241       returnedSupportDescriptionChar = new char[lenDescription];
2242
2243       returnedSupportNameChar = strcpy(returnedSupportNameChar,"Copy of ");
2244       returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
2245       returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Copy of ");
2246       returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2247                                               (Supports[0]->getDescription()).c_str());
2248
2249       returnedSupportName = string(returnedSupportNameChar);
2250       returnedSupportDescription = string(returnedSupportDescriptionChar);
2251
2252       returnedSupport->setName(returnedSupportName);
2253       returnedSupport->setDescription(returnedSupportDescription);
2254
2255       delete [] returnedSupportNameChar;
2256       delete [] returnedSupportDescriptionChar;
2257     }
2258   else
2259     {
2260       SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2261       returnedSupport = new SUPPORT(*obj);
2262
2263       int lenName = strlen((Supports[0]->getName()).c_str()) + 16 + 1;
2264       int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 16 + 1;
2265
2266       for(int i = 1;i<size;i++)
2267         {
2268           obj = const_cast <SUPPORT *> (Supports[i]);
2269           returnedSupport->intersecting(obj);
2270
2271           if (i == (size-1))
2272             {
2273               lenName = lenName + 5 + strlen((Supports[i]->getName()).c_str());
2274               lenDescription = lenDescription + 5 +
2275                 strlen((Supports[i]->getDescription()).c_str());
2276             }
2277           else
2278             {
2279               lenName = lenName + 2 + strlen((Supports[i]->getName()).c_str());
2280               lenDescription = lenDescription + 2 +
2281                 strlen((Supports[i]->getDescription()).c_str());
2282             }
2283         }
2284       if(returnedSupport != (SUPPORT *) NULL)
2285         {
2286           returnedSupportNameChar = new char[lenName];
2287           returnedSupportDescriptionChar = new char[lenDescription];
2288
2289           returnedSupportNameChar = strcpy(returnedSupportNameChar,
2290                                            "Intersection of ");
2291           returnedSupportDescriptionChar =
2292             strcpy(returnedSupportDescriptionChar,"Intersection of ");
2293
2294           returnedSupportNameChar = strcat(returnedSupportNameChar,
2295                                            (Supports[0]->getName()).c_str());
2296           returnedSupportDescriptionChar =
2297             strcat(returnedSupportDescriptionChar,
2298                    (Supports[0]->getDescription()).c_str());
2299
2300           for(int i = 1;i<size;i++)
2301             {
2302               if (i == (size-1))
2303                 {
2304                   returnedSupportNameChar = strcat(returnedSupportNameChar,
2305                                                    " and ");
2306                   returnedSupportDescriptionChar =
2307                     strcat(returnedSupportDescriptionChar," and ");
2308
2309                   returnedSupportNameChar =
2310                     strcat(returnedSupportNameChar,
2311                            (Supports[i]->getName()).c_str());
2312                   returnedSupportDescriptionChar =
2313                     strcat(returnedSupportDescriptionChar,
2314                            (Supports[i]->getDescription()).c_str());
2315                 }
2316               else
2317                 {
2318                   returnedSupportNameChar = strcat(returnedSupportNameChar,
2319                                                    ", ");
2320                   returnedSupportNameChar =
2321                     strcat(returnedSupportNameChar,
2322                            (Supports[i]->getName()).c_str());
2323
2324                   returnedSupportDescriptionChar =
2325                     strcat(returnedSupportDescriptionChar,", ");
2326                   returnedSupportDescriptionChar =
2327                     strcat(returnedSupportDescriptionChar,
2328                            (Supports[i]->getDescription()).c_str());
2329                 }
2330             }
2331
2332           returnedSupportName = string(returnedSupportNameChar);
2333           returnedSupport->setName(returnedSupportName);
2334
2335           returnedSupportDescription = string(returnedSupportDescriptionChar);
2336           returnedSupport->setDescription(returnedSupportDescription);
2337
2338           delete [] returnedSupportNameChar;
2339           delete [] returnedSupportDescriptionChar;
2340         }
2341     }
2342
2343   END_OF_MED(LOC);
2344   return returnedSupport;
2345 }
2346 /*!
2347 @}
2348 */
2349
2350 // internal helper type
2351 struct _cell
2352 {
2353     std::vector<int> groups;
2354     MED_EN::medGeometryElement geometricType;
2355 };
2356
2357 /*!\addtogroup MESH_families
2358 @{
2359 */
2360
2361 /*! Retrieves the group named \a name.
2362 The method browses all the entities in order to find the group.
2363 If two groups with the same name coexist, the first one found will be
2364 returned. If no group with the correct name is found, the method throws
2365 an exception.
2366  */
2367 const GROUP* MESH::getGroup(const string& name) const  throw (MEDEXCEPTION)
2368 {
2369         const vector<GROUP*>* group_vectors [4]={&_groupNode, &_groupEdge,&_groupFace,&_groupCell};
2370         for (int ientity=0;ientity<4;ientity++)
2371                 for (int igroup=0; igroup< group_vectors[ientity]->size();igroup++)
2372                         {
2373                                 const vector<GROUP*>& group_vect = *group_vectors[ientity];
2374                                 GROUP* group=group_vect[igroup];
2375                                 if (group->getName()==name)
2376                                         return group;
2377                         }
2378         cerr << "MESH::getGroup("<<name<<") : group "<<name <<" was not found"<<endl;
2379         throw MEDEXCEPTION("MESH::getGroup(name) : name not found");
2380 }
2381 /*!
2382 @}
2383 */
2384
2385 const GROUP* MESH::getGroup(MED_EN::medEntityMesh entity, int i) const
2386 {
2387   const char * LOC = "MESH::getGroup(medEntityMesh entity, int i) : ";
2388   if (i<=0)
2389     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"argument i must be > 0"));
2390   vector<GROUP*> Group;
2391   switch (entity) {
2392   case MED_EN::MED_NODE : {
2393     Group = _groupNode;
2394     break;
2395   }
2396   case MED_EN::MED_CELL : {
2397     Group = _groupCell;
2398     break;
2399   }
2400   case MED_EN::MED_FACE : {
2401     Group = _groupFace;
2402     break;
2403   }
2404   case MED_EN::MED_EDGE : {
2405     Group = _groupEdge;
2406     break;
2407   }
2408   default :
2409     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Unknown entity"));
2410   }
2411   if (i>(int)Group.size())
2412     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"argument i="<<i<<" must be <= _numberOfGroups="<<Group.size()));
2413   return Group[i-1];
2414 }
2415
2416
2417 /*!
2418 \addtogroup MESH_families
2419 @{
2420 */
2421 /*! Returns the groups of type \a entity present in the mesh as a vector of pointers. The GROUP class inheriting from the SUPPORT class, the
2422 methods that can be used on these groups are explained in the related section.*/
2423 const vector<GROUP*> MESH::getGroups(MED_EN::medEntityMesh entity) const
2424 {
2425   switch (entity) {
2426   case MED_EN::MED_NODE :
2427     return _groupNode;
2428   case MED_EN::MED_CELL :
2429     return _groupCell;
2430   case MED_EN::MED_FACE :
2431     return _groupFace;
2432   case MED_EN::MED_EDGE :
2433     return _groupEdge;
2434   default :
2435     throw MEDEXCEPTION("MESH::getGroups : Unknown entity");
2436   }
2437 }
2438 /*!
2439 @}
2440 */
2441
2442 /*!
2443   Create groups from families.
2444
2445   It is used to create groups that have only one family
2446   for meshes that come from codes that use families instead
2447   of groups to define a subregion.
2448 */
2449 void MESH::createGroups()
2450 {
2451   for (medEntityMesh entity=MED_CELL; entity!=MED_ALL_ENTITIES; ++entity)
2452     {
2453       // make myFamilies points to the member corresponding to entity
2454       vector<FAMILY*>* myFamilies;
2455       vector<GROUP*>* myGroups;
2456       switch ( entity )
2457         {
2458         case MED_CELL :
2459           myFamilies = & _familyCell;
2460           myGroups = & _groupCell;
2461           break;
2462         case MED_FACE :
2463           myFamilies = & _familyFace;
2464           myGroups = & _groupFace;
2465           break;
2466         case MED_EDGE :
2467           myFamilies = & _familyEdge;
2468           myGroups = & _groupEdge;
2469           break;
2470         case MED_NODE :
2471           myFamilies = & _familyNode;
2472           myGroups = & _groupNode;
2473           break;
2474         }
2475
2476
2477       for (int i=0; i< myFamilies->size(); i++)
2478         {
2479           list <FAMILY*> fam_list;
2480           fam_list.push_back((*myFamilies)[i]);
2481           //creates a group with the family name and only one family
2482           GROUP* group=new GROUP((*myFamilies)[i]->getName(),fam_list);
2483           (*myGroups).push_back(group);
2484         }
2485     }
2486 }
2487 // Create families from groups
2488 void MESH::createFamilies()
2489 {
2490     int idFamNode = 0; // identifier for node families
2491     int idFamElement = 0; // identifier for cell, face or edge families
2492
2493     // Main loop on mesh's entities
2494     for (medEntityMesh entity=MED_CELL; entity!=MED_ALL_ENTITIES; ++entity)
2495     {
2496         int numberofgroups = getNumberOfGroups(entity);
2497         if(!numberofgroups)
2498             continue; // no groups for this entity
2499
2500         vector< vector<FAMILY*> > whichFamilyInGroup(numberofgroups); // this container is used to update groups at the end
2501
2502         // make myFamilies points to the member corresponding to entity
2503         vector<FAMILY*>* myFamilies;
2504         switch ( entity )
2505         {
2506             case MED_CELL :
2507                 myFamilies = & _familyCell;
2508                 break;
2509             case MED_FACE :
2510                 myFamilies = & _familyFace;
2511                 break;
2512             case MED_EDGE :
2513                 myFamilies = & _familyEdge;
2514                 break;
2515             case MED_NODE :
2516                 myFamilies = & _familyNode;
2517                 break;
2518         }
2519
2520         vector<GROUP*> myGroups=getGroups(entity); // get a copy of the groups ptr for the entity
2521         // get a copy of the (old) family ptrs before clearing
2522         vector<FAMILY*> myOldFamilies=getFamilies(entity);
2523         myFamilies->clear();
2524
2525
2526         // 1 - Create a vector containing for each cell (of the entity) an information structure
2527         //     giving geometric type and the groups it belong to
2528
2529         med_int numberOfTypes=getNumberOfTypesWithPoly(entity);
2530         medGeometryElement* geometricTypes=_connectivity->getGeometricTypesWithPoly(entity); // pb avec entity=MED_NODE???
2531         med_int numberOfCells=getNumberOfElementsWithPoly(entity, MED_ALL_ELEMENTS);  // total number of cells for that entity
2532         SCRUTE_MED(numberOfTypes);
2533         SCRUTE_MED(numberOfCells);
2534         vector< _cell > tab_cell(numberOfCells);
2535         vector< _cell >::iterator cell = tab_cell.begin();
2536         for(med_int t=0; t!=numberOfTypes; ++t) {
2537           int nbCellsOfType = getNumberOfElementsWithPoly(entity,geometricTypes[t]);
2538           for(int n=0; n!=nbCellsOfType; ++n, ++cell)
2539             cell->geometricType=geometricTypes[t];
2540         }
2541         delete [] geometricTypes;
2542
2543         // 2 - Scan cells in groups and update in tab_cell the container of groups a cell belong to
2544
2545         for (unsigned g=0; g!=myGroups.size(); ++g)
2546         {
2547             // scan cells that belongs to the group
2548             const int* groupCells=myGroups[g]->getnumber()->getValue();
2549             int nbCells=myGroups[g]->getnumber()->getLength();
2550             for(int c=0; c!=nbCells; ++c)
2551                 tab_cell[groupCells[c]-1].groups.push_back(g);
2552         }
2553
2554
2555         // 3 - Scan the cells vector, genarate family name, and create a map associating the family names
2556         //     whith the vector of contained cells
2557
2558         map< string,vector<int> > tab_families;
2559         map< string,vector<int> >::iterator fam;
2560         for(int n=0; n!=numberOfCells; ++n)
2561         {
2562             ostringstream key; // to generate the name of the family
2563             key << "FAM";
2564             if(tab_cell[n].groups.empty()) // this cell don't belong to any group
2565                 key << "_NONE" << entity;
2566
2567             for(vector<int>::const_iterator it=tab_cell[n].groups.begin(); it!=tab_cell[n].groups.end(); ++it)
2568             {
2569                 string groupName=myGroups[*it]->getName();
2570                 if(groupName.empty())
2571                     key << "_G" << *it;
2572                 else
2573                     key << "_" << groupName;
2574             }
2575
2576             tab_families[key.str()].push_back(n+1); // fill the vector of contained cells associated whith the family
2577         }
2578
2579
2580         // 4 - Scan the family map, create MED Families, check if it already exist.
2581
2582         for( fam=tab_families.begin(); fam!=tab_families.end(); ++fam)
2583         {
2584             vector<medGeometryElement> tab_types_geometriques;
2585             medGeometryElement geometrictype=MED_NONE;
2586             vector<int> tab_index_types_geometriques;
2587             vector<int> tab_nombres_elements;
2588             if ( fam->second.empty() )
2589               continue; // it is just a truncated long family name
2590
2591             // scan family cells and fill the tab that are needed by the create a MED FAMILY
2592             for( int i=0; i!=fam->second.size(); ++i)
2593             {
2594                 int ncell=fam->second[i]-1;
2595                 if(tab_cell[ncell].geometricType != geometrictype)
2596                 {
2597                     // new geometric type -> we store it and complete index tabs
2598                     if(!tab_index_types_geometriques.empty())
2599                         tab_nombres_elements.push_back(i+1-tab_index_types_geometriques.back());
2600                     tab_types_geometriques.push_back( (geometrictype=tab_cell[ncell].geometricType));
2601                     tab_index_types_geometriques.push_back(i+1);
2602                 }
2603             }
2604             // store and complete index tabs for the last geometric type
2605             tab_nombres_elements.push_back(fam->second.size()+1-tab_index_types_geometriques.back());
2606             tab_index_types_geometriques.push_back(fam->second.size()+1);
2607
2608             // family name sould not be longer than MED_TAILLE_NOM
2609             string famName = fam->first;
2610             if ( famName.size() > MED_TAILLE_NOM ) {
2611               // try to cut off "FAM_" from the head
2612               if ( famName.size() - 4 <= MED_TAILLE_NOM ) {
2613                 famName = famName.substr(4);
2614               }
2615               else { // try to make a unique name by cutting off char by char from the tail
2616                 famName = famName.substr(0, MED_TAILLE_NOM);
2617                 map< string,vector<int> >::iterator foundName = tab_families.find( famName );
2618                 while ( !famName.empty() &&
2619                         ( foundName != tab_families.end() || famName[ famName.size()-1 ] == ' ' ))
2620                 {
2621                   famName = famName.substr( 0, famName.size() - 1 );
2622                   foundName = tab_families.find( famName );
2623                 }
2624               }
2625               tab_families[ famName ]; // add a new name in the table to assure uniqueness
2626             }
2627
2628             // create an empty MED FAMILY and fill it with the tabs we constructed
2629             FAMILY* newFam = new FAMILY();
2630             //newFam->setTotalNumberOfElements(fam->second.size());
2631             newFam->setName(famName);
2632             newFam->setMeshDirectly(this);
2633             newFam->setNumberOfGeometricType(tab_types_geometriques.size());
2634             newFam->setGeometricType(&tab_types_geometriques[0]); // we know the tab is not empy
2635             newFam->setNumberOfElements(&tab_nombres_elements[0]);
2636             newFam->setNumber(&tab_index_types_geometriques[0],&fam->second[0]);
2637             newFam->setEntity(entity);
2638             newFam->setAll(false);
2639
2640             int idFam = 0;
2641
2642             switch ( entity )
2643               {
2644               case MED_NODE :
2645                 ++idFamNode;
2646                 idFam = idFamNode;
2647                 break;
2648               case MED_CELL:
2649                 ++idFamElement;
2650                 idFam = -idFamElement;
2651                 break;
2652               case MED_FACE :
2653                 ++idFamElement;
2654                 idFam = -idFamElement;
2655                 break;
2656               case MED_EDGE :
2657                 ++idFamElement;
2658                 idFam = -idFamElement;
2659                 break;
2660               }
2661
2662             newFam->setIdentifier(idFam);
2663
2664             // Update links between families and groups
2665
2666             int ncell1=fam->second[0]-1;  // number of first cell in family
2667             int numberOfGroups=tab_cell[ncell1].groups.size(); // number of groups in the family
2668             if(numberOfGroups)
2669             {
2670                 newFam->setNumberOfGroups(numberOfGroups);
2671                 string * groupNames=new string[numberOfGroups];
2672
2673                 // iterate on the groups the family belongs to
2674                 vector<int>::const_iterator it=tab_cell[ncell1].groups.begin();
2675                 for(int ng=0 ; it!=tab_cell[ncell1].groups.end(); ++it, ++ng)
2676                 {
2677                     whichFamilyInGroup[*it].push_back(newFam);
2678                     groupNames[ng]=myGroups[*it]->getName();
2679                 }
2680                 newFam->setGroupsNames(groupNames);
2681             }
2682
2683             MESSAGE_MED("  MESH::createFamilies() entity " << entity <<
2684                         " size " << myFamilies->size());
2685
2686             myFamilies->push_back(newFam);
2687         }
2688
2689         // delete old families
2690         for (unsigned int i=0;i<myOldFamilies.size();i++)
2691             delete myOldFamilies[i] ;
2692
2693         // update references in groups
2694         for (unsigned int i=0;i<myGroups.size();i++)
2695         {
2696             myGroups[i]->setNumberOfFamilies(whichFamilyInGroup[i].size());
2697             myGroups[i]->setFamilies(whichFamilyInGroup[i]);
2698         }
2699
2700     // re-scan the cells vector, and fill the family vector with cells.
2701     // creation of support, check if it already exist.
2702     }
2703 }
2704
2705 int MESH::getElementContainingPoint(const double *coord)
2706 {
2707   if(_spaceDimension==3)
2708     {
2709       Meta_Wrapper<3> fromWrapper(getNumberOfNodes(),const_cast<double *>(getCoordinates(MED_FULL_INTERLACE)),
2710                                                         const_cast<CONNECTIVITY *>(getConnectivityptr()));
2711       Meta_Wrapper<3> toWrapper(1,const_cast<double *>(coord));
2712       Meta_Mapping<3> mapping(&fromWrapper,&toWrapper);
2713       mapping.Cree_Mapping(1);
2714       return mapping.Get_Mapping()[0]+1;
2715     }
2716   else if(_spaceDimension==2)
2717     {
2718       Meta_Wrapper<2> fromWrapper(getNumberOfNodes(),const_cast<double *>(getCoordinates(MED_FULL_INTERLACE)),
2719                                                         const_cast<CONNECTIVITY *>(getConnectivityptr()));
2720       Meta_Wrapper<2> toWrapper(1,const_cast<double *>(coord));
2721       Meta_Mapping<2> mapping(&fromWrapper,&toWrapper);
2722       mapping.Cree_Mapping(1);
2723       return mapping.Get_Mapping()[0]+1;
2724       }
2725   else
2726     throw MEDEXCEPTION("MESH::getElementContainingPoint : invalid _spaceDimension must be equal to 2 or 3 !!!");
2727 }
2728
2729 //! Converts MED_CELL connectivity to polyhedra connectivity
2730 //! Converts MED_FACE connectivity to polygon connectivity
2731 //! Wil work only for 3D meshes with nodal connectivity
2732
2733 void MESH::convertToPoly()
2734 {
2735   if (getMeshDimension()!=3) return;
2736
2737   CONNECTIVITY* newpolygonconnectivity = new CONNECTIVITY(MED_EN::MED_FACE);
2738   CONNECTIVITY* newpolyhedraconnectivity = new CONNECTIVITY(MED_EN::MED_CELL);
2739
2740   {
2741     ////////////////////////////////////////////:
2742     // First step : Treating polygons connectivity
2743     ///////////////////////////////////////////
2744
2745     const int* oldconn = getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL, MED_EN::MED_FACE, MED_EN::MED_ALL_ELEMENTS);
2746
2747     const int* oldconnindex= getConnectivityIndex(MED_EN::MED_NODAL,MED_EN::MED_FACE);
2748     int oldnbface = getNumberOfElements(MED_EN::MED_FACE,MED_EN::MED_ALL_ELEMENTS);
2749     const int* oldconnpoly =0;
2750     const int* oldconnpolyindex =0;
2751     if(existPolygonsConnectivity(MED_EN::MED_NODAL, MED_EN::MED_FACE))
2752     {
2753      oldconnpoly = getPolygonsConnectivity(MED_EN::MED_NODAL, MED_EN::MED_FACE);
2754      oldconnpolyindex = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,MED_EN::MED_FACE);
2755      }
2756     int oldnbtotalface = getNumberOfElementsWithPoly(MED_EN::MED_FACE,MED_EN::MED_ALL_ELEMENTS);
2757
2758     int nbnodes=0;
2759
2760     if (oldconnindex !=0)
2761       nbnodes += oldconnindex[oldnbface] -1 ;
2762
2763     if (oldconnpolyindex !=0)
2764       nbnodes+= oldconnpolyindex[oldnbtotalface-oldnbface]-1;
2765
2766     int* newconn = new int[nbnodes];
2767     int* newconnindex= new int [oldnbtotalface+1];
2768
2769     //copying classical types connectivity
2770     memcpy(newconn, oldconn, sizeof(int)*(oldconnindex[oldnbface]-1) );
2771
2772     //copying poly types connectivity
2773     if(oldconnpoly)
2774       memcpy (newconn+oldconnindex[oldnbface]-1, oldconnpoly, sizeof(int)*(oldconnpolyindex[oldnbtotalface-oldnbface]-1) );
2775
2776     newconnindex[0]=1;
2777     for (int i=0; i<oldnbface;i++)
2778       newconnindex[i+1]=newconnindex[i]+oldconnindex[i+1]-oldconnindex[i];
2779     for (int i=oldnbface; i<oldnbtotalface;i++)
2780       newconnindex[i+1]=newconnindex[i]+
2781         oldconnpolyindex[i-oldnbface+1]-oldconnpolyindex[i-oldnbface];
2782
2783
2784     newpolygonconnectivity->setPolygonsConnectivity(MED_EN::MED_NODAL,
2785                                              MED_EN::MED_FACE,
2786                                              newconn,
2787                                              newconnindex,
2788                                              nbnodes,
2789                                              oldnbtotalface);
2790     //    _connectivity->setConstituent(newconnectivity);
2791   }
2792   ///////////////////////////////////////////
2793   // 2nd step : Treating polyhedra connectivity
2794   //////////////////////////////////////////
2795   {
2796
2797     const int* oldconn = getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL, MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS);
2798
2799     const int* oldconnindex= getConnectivityIndex(MED_EN::MED_NODAL,MED_EN::MED_CELL);
2800     int oldnbelem = getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
2801     const int* oldconnpoly = 0;
2802     const int* oldconnpolyindex = 0;
2803     const int* oldfaceindex = 0;
2804     if(existPolyhedronConnectivity(MED_EN::MED_NODAL, MED_EN::MED_CELL))
2805     {
2806        oldconnpoly = getPolyhedronConnectivity(MED_EN::MED_NODAL);
2807        oldconnpolyindex = getPolyhedronIndex(MED_EN::MED_NODAL);
2808        oldfaceindex =  getPolyhedronFacesIndex();
2809     }
2810     const MED_EN::medGeometryElement* oldtypes = getTypes(MED_EN::MED_CELL);
2811     int nboldtypes=getNumberOfTypes(MED_EN::MED_CELL);
2812     int nboldpolyhedra=getNumberOfPolyhedron();
2813     int oldnbtotalelem = getNumberOfElementsWithPoly(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
2814
2815     int nbnodes=0;
2816
2817     if (oldconnindex !=0)
2818       nbnodes += oldconnindex[oldnbelem] -1 ;
2819
2820     if (oldconnpolyindex !=0)
2821       nbnodes+= oldconnpolyindex[oldnbtotalelem-oldnbelem]-1;
2822
2823     //computing number of faces
2824     int nbfaces=0;
2825     //    first part : number of faces for the classical types
2826     for (int itype=0; itype<nboldtypes; itype++)
2827       {
2828         MED_EN::medGeometryElement type = oldtypes[itype];
2829         MEDMEM::CELLMODEL cellmodel(type);
2830         int nb_elems=getNumberOfElements(MED_EN::MED_CELL,type);
2831         int nbfacespertype = cellmodel.getNumberOfConstituents(1);
2832         nbfaces+=nb_elems*nbfacespertype;
2833       }
2834     //   second part : number of faces for the polyhedra
2835     nbfaces += getNumberOfPolyhedronFaces();
2836
2837     //allocating tables for new connectivity
2838   vector<int> newconn;
2839   vector<int> newconnindex(1,1);
2840   vector<int> newfaceindex(1,1);
2841
2842   for (int itype=0; itype<nboldtypes; itype++)
2843     {
2844       MED_EN::medGeometryElement type = oldtypes[itype];
2845       MEDMEM::CELLMODEL cellmodel(type);
2846       int nb_elems=getNumberOfElements(MED_EN::MED_CELL,type);
2847       int nbfacespertype = cellmodel.getNumberOfConstituents(1);
2848         for (int ielem = 0; ielem<nb_elems; ielem++)
2849           {
2850             for (int iface =0; iface< nbfacespertype; iface ++)
2851               {
2852                 //local conn contains the local nodal connectivity for the iface-th face of type type
2853                 const int* local_conn = cellmodel.getNodesConstituent(1,iface+1); // iface+1 for MED numbering
2854                 MED_EN::medGeometryElement facetype = cellmodel.getConstituentType(1,iface+1);
2855                 int nbface_nodes=facetype%100;
2856                 for ( int inode=0; inode<nbface_nodes;inode++)
2857                   {
2858                     newconn.push_back(oldconn[oldconnindex[newconnindex.size()-1]-1+local_conn[inode]-1]);
2859                   }
2860                 newfaceindex.push_back(newfaceindex[newfaceindex.size()-1]+nbface_nodes);
2861               }
2862             newconnindex.push_back(newconnindex[newconnindex.size()-1]+nbfacespertype);
2863           }
2864       }
2865
2866   for (int i=0; i<nboldpolyhedra; i++)
2867     {
2868       newconnindex.push_back(newconnindex[newconnindex.size()-1]+oldconnpolyindex[i+1]-oldconnpolyindex[i]);
2869     }
2870   if(oldconnpolyindex)
2871   {
2872     for (int i=0; i<oldconnpolyindex[nboldpolyhedra]-1;i++)
2873       {
2874         newfaceindex.push_back(newfaceindex[newfaceindex.size()-1]+oldfaceindex[i+1]-oldfaceindex[i]);
2875       }
2876     for (int i=0; i< oldfaceindex[oldconnpolyindex[nboldpolyhedra]-1]-1; i++)
2877       newconn.push_back(oldconnpoly[i]);
2878   }
2879   //  memcpy(newconn_ptr,oldconnpoly,sizeof(int)*(oldfaceindex[oldconnpoly[nboldpolyhedra]-1]-1));
2880
2881
2882   newpolyhedraconnectivity->setPolyhedronConnectivity(MED_EN::MED_NODAL,
2883                                              &newconn[0],
2884                                              &newconnindex[0],
2885                                              newfaceindex[newfaceindex.size()-1]-1,
2886                                              newconnindex.size()-1,
2887                                              &newfaceindex[0],
2888                                              newfaceindex.size()-1);
2889
2890   newpolyhedraconnectivity->setEntityDimension(3);
2891
2892   delete _connectivity;
2893
2894  _connectivity=newpolyhedraconnectivity;
2895  _connectivity->setConstituent(newpolygonconnectivity);
2896
2897   }
2898 }
2899
2900 vector< vector<double> > MESH::getBoundingBox() const
2901 {
2902   const double *myCoords=_coordinate->getCoordinates(MED_EN::MED_FULL_INTERLACE);
2903   vector< vector<double> > ret(2);
2904   int i,j;
2905   ret[0].resize(_spaceDimension);
2906   ret[1].resize(_spaceDimension);
2907   for(i=0;i<_spaceDimension;i++)
2908     {
2909       ret[0][i]=1.e300;
2910       ret[1][i]=-1.e300;
2911     }
2912   for(i=0;i<_coordinate->getNumberOfNodes();i++)
2913     {
2914       for(j=0;j<_spaceDimension;j++)
2915         {
2916           double tmp=myCoords[i*_spaceDimension+j];
2917           if(tmp<ret[0][j])
2918             ret[0][j]=tmp;
2919           if(tmp>ret[1][j])
2920             ret[1][j]=tmp;
2921         }
2922     }
2923   return ret;
2924 }
2925
2926 //Presently disconnected in C++
2927 void MESH::addReference() const
2928 {
2929 }
2930
2931 //Presently disconnected in C++
2932 void MESH::removeReference() const
2933 {
2934 }