src/MEDMEM/MEDMEM_Mesh.cxx

   1 /*
   2  File Mesh.cxx
   3  $Header$
   4 */
   5
   6 #include <math.h>
   7
   8 #include <list>
   9 #include <map>
  10
  11 #include "MEDMEM_Field.hxx"
  12 #include "MEDMEM_Mesh.hxx"
  13
  14 #include "MEDMEM_Family.hxx"
  15 #include "MEDMEM_Group.hxx"
  16 #include "MEDMEM_Connectivity.hxx"
  17 #include "MEDMEM_CellModel.hxx"
  18
  19 //update Families with content list
  20 //int family_count(int family_number, int count, int * entities_number, int * entities_list) ;
  21
  22 // ------- Drivers Management Part
  23
  24 // MESH::INSTANCE_DE<MED_MESH_RDONLY_DRIVER> MESH::inst_med_rdonly ;
  25 //const MESH::INSTANCE * const MESH::instances[] = { &MESH::inst_med_rdonly , &MESH::inst_med_rdwr } ;
  26
  27 MESH::INSTANCE_DE<MED_MESH_RDWR_DRIVER>   MESH::inst_med   ;
  28 const MESH::INSTANCE * const MESH::instances[] =   {  &MESH::inst_med } ;
  29
  30 /*! Add a MESH driver of type (MED_DRIVER, ....) associated with file <fileName>. The meshname used in the file
  31     is  <driverName>. addDriver returns an int handler. */
  32 int MESH::addDriver(driverTypes driverType,
  33                     const string & fileName="Default File Name.med",const string & driverName="Default Mesh Name") {
  34
  35   const char * LOC = "MESH::addDriver(driverTypes driverType, const string & fileName=\"Default File Name.med\",const string & driverName=\"Default Mesh Name\") : ";
  36
  37   GENDRIVER * driver;
  38   int current;
  39
  40   BEGIN_OF(LOC);
  41
  42   driver = instances[driverType]->run(fileName, this) ;
  43   _drivers.push_back(driver);
  44   current = _drivers.size()-1;
  45
  46   _drivers[current]->setMeshName(driverName);
  47   return current;
  48
  49   END_OF(LOC);
  50
  51 }
  52
  53 /*! Add an existing MESH driver. */
  54 int  MESH::addDriver(GENDRIVER & driver) {
  55   const char * LOC = "MESH::addDriver(GENDRIVER &) : ";
  56   BEGIN_OF(LOC);
  57
  58   // A FAIRE VERIFIER QUE LE DRIVER EST DE TYPE MESH !!
  59   _drivers.push_back(&driver);
  60   return _drivers.size()-1;
  61
  62   END_OF(LOC);
  63 }
  64
  65 /*! Remove an existing MESH driver. */
  66 void MESH::rmDriver (int index=0) {
  67   const char * LOC = "MESH::rmDriver (int index=0): ";
  68   BEGIN_OF(LOC);
  69
  70   if ( _drivers[index] ) {
  71     //_drivers.erase(&_drivers[index]);
  72     // why not ????
  73     MESSAGE ("detruire");
  74   }
  75   else
  76     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
  77                                      << "The index given is invalid, index must be between  0 and  |"
  78                                      << _drivers.size()
  79                                      )
  80                           );
  81
  82   END_OF(LOC);
  83
  84 };
  85
  86 // ------ End of Drivers Management Part
  87
  88
  89 void MESH::init() {
  90
  91   string        _name = "NOT DEFINED"; // A POSITIONNER EN FCT DES IOS ?
  92
  93   _numberOfMEDNodeFamily =          MED_INVALID;
  94   _MEDArrayNodeFamily    = (int * ) NULL; // SOLUTION TEMPORAIRE
  95   _numberOfMEDCellFamily = (int * ) NULL;
  96   _numberOfMEDFaceFamily = (int * ) NULL;
  97   _numberOfMEDEdgeFamily = (int * ) NULL;
  98   _MEDArrayCellFamily    = (int **) NULL; // SOLUTION TEMPORAIRE
  99   _MEDArrayFaceFamily    = (int **) NULL; // SOLUTION TEMPORAIRE
 100   _MEDArrayEdgeFamily    = (int **) NULL; // SOLUTION TEMPORAIRE
 101
 102   COORDINATE   * _coordinate   = (COORDINATE   *) NULL;
 103   CONNECTIVITY * _connectivity = (CONNECTIVITY *) NULL;
 104
 105   _spaceDimension        =          MED_INVALID; // 0 ?!?
 106   _meshDimension         =          MED_INVALID;
 107   _numberOfNodes         =          MED_INVALID;
 108
 109   _numberOfNodesFamilies =          0; // MED_INVALID ?!?
 110   _numberOfCellsFamilies =          0;
 111   _numberOfFacesFamilies =          0;
 112   _numberOfEdgesFamilies =          0;
 113
 114   _numberOfNodesGroups   =          0;  // MED_INVALID ?!?
 115   _numberOfCellsGroups   =          0;
 116   _numberOfFacesGroups   =          0;
 117   _numberOfEdgesGroups   =          0;
 118
 119 };
 120
 121 /*! Create an empty MESH. */
 122 MESH::MESH() {
 123   init();
 124 };
 125
 126 MESH::MESH(const MESH &m)
 127 {
 128   // VERIFIER QUE TS LES OPERATEURS DE RECOPIE DES ATTRIBUTS
 129   // SONT CORRECTS.
 130   *this = m;
 131 }
 132
 133 MESH::~MESH() {
 134
 135   if ( _MEDArrayNodeFamily    != (int * ) NULL) delete [] _MEDArrayNodeFamily; // SOLUTION TEMPORAIRE
 136   if ( _numberOfMEDCellFamily != (int * ) NULL) delete [] _numberOfMEDCellFamily;
 137   if ( _numberOfMEDFaceFamily != (int * ) NULL) delete [] _numberOfMEDFaceFamily;
 138   if ( _numberOfMEDEdgeFamily != (int * ) NULL) delete [] _numberOfMEDEdgeFamily;
 139   // IL FAUT FAIRE UNE BOUCLE DE DESALLOCATION
 140   if ( _MEDArrayCellFamily    != (int **) NULL) delete [] _MEDArrayCellFamily; // SOLUTION TEMPORAIRE
 141   if ( _MEDArrayFaceFamily    != (int **) NULL) delete [] _MEDArrayFaceFamily; // SOLUTION TEMPORAIRE
 142   if ( _MEDArrayEdgeFamily    != (int **) NULL) delete [] _MEDArrayEdgeFamily; // SOLUTION TEMPORAIRE
 143
 144   if (_coordinate != NULL) delete _coordinate ;
 145   if (_connectivity != NULL) delete _connectivity ;
 146   int size ;
 147   size = _familyNode.size() ;
 148   for (int i=0;i<size;i++)
 149     delete _familyNode[i] ;
 150   size = _familyCell.size() ;
 151   for (int i=0;i<size;i++)
 152     delete _familyCell[i] ;
 153   size = _familyFace.size() ;
 154   for (int i=0;i<size;i++)
 155     delete _familyFace[i] ;
 156   size = _familyEdge.size() ;
 157   for (int i=0;i<size;i++)
 158     delete _familyEdge[i] ;
 159   size = _groupNode.size() ;
 160   for (int i=0;i<size;i++)
 161     delete _groupNode[i] ;
 162   size = _groupCell.size() ;
 163   for (int i=0;i<size;i++)
 164     delete _groupCell[i] ;
 165   size = _groupFace.size() ;
 166   for (int i=0;i<size;i++)
 167     delete _groupFace[i] ;
 168   size = _groupEdge.size() ;
 169   for (int i=0;i<size;i++)
 170     delete _groupEdge[i] ;
 171
 172 }
 173
 174 MESH & MESH::operator=(const MESH &m)
 175 {
 176
 177   //  A FAIRE.........
 178
 179   // ATTENTION CET OPERATEUR DE RECOPIE EST DANGEREUX POUR LES
 180   // POINTEURS : ex : nodal_connectivity ???? EXPRES ????
 181 //    _drivers = m._drivers;
 182
 183 //    space_dimension=m.space_dimension;
 184 //    mesh_dimension=m.mesh_dimension;
 185
 186 //    nodes_count=m.nodes_count;
 187 //    coordinates=m.coordinates;
 188 //    coordinates_name=m.coordinates_name ;
 189 //    coordinates_unit=m.coordinates_unit ;
 190 //    nodes_numbers=m.nodes_numbers ;
 191
 192 //    cells_types_count=m.cells_types_count;
 193 //    cells_type=m.cells_type;
 194 //    cells_count=m.cells_count;
 195 //    nodal_connectivity=m.nodal_connectivity;
 196
 197 //    nodes_families_count=m.nodes_families_count;
 198 //    nodes_Families=m.nodes_Families;
 199
 200 //    cells_families_count=m.cells_families_count;
 201 //    cells_Families=m.cells_Families;
 202
 203 //    maximum_cell_number_by_node = m.maximum_cell_number_by_node;
 204 //    if (maximum_cell_number_by_node > 0)
 205 //      {
 206 //        reverse_nodal_connectivity = m.reverse_nodal_connectivity;
 207 //        reverse_nodal_connectivity_index = m.reverse_nodal_connectivity_index ;
 208 //      }
 209   return *this;
 210 }
 211
 212 /*! Create a MESH object using a MESH driver of type (MED_DRIVER, ....) associated with file <fileName>.
 213   The meshname <driverName> must already exists in the file.*/
 214 MESH::MESH(driverTypes driverType, const string &  fileName="", const string &  driverName="") {
 215   const char * LOC ="MESH::MESH(driverTypes driverType, const string &  fileName="", const string &  driverName="") : ";
 216
 217   int current;
 218
 219   BEGIN_OF(LOC);
 220
 221   init();
 222
 223   current = addDriver(driverType,fileName,driverName);
 224   switch(_drivers[current]->getAccessMode() ) {
 225   case MED_RDONLY : {
 226     MESSAGE("MESH::MESH(driverTypes driverType, .....) : driverType must have a MED_RDWR accessMode");
 227     rmDriver(current);
 228     break;}
 229   default : {
 230   }
 231   }
 232   _drivers[current]->open();
 233   _drivers[current]->read();
 234   _drivers[current]->close();
 235   END_OF(LOC);
 236 };
 237
 238
 239 //  Node MESH::Donne_Barycentre(const Maille &m) const
 240 //  {
 241 //    Node N=node[m[0]];
 242 //    for (int i=1;i<m.donne_nbr_sommets();i++) N=N+node[m[i]];
 243 //    N=N*(1.0/m.donne_nbr_sommets());
 244 //    return N;
 245 //  }
 246
 247 ostream & operator<<(ostream &os, MESH &myMesh)
 248 {
 249 //    int space_dimension = myMesh.get_space_dimension();
 250 //    os << "Space dimension "<< space_dimension << endl ;
 251 //    int nodes_count = myMesh.get_nodes_count() ;
 252 //    os << "Nodes count : "<< nodes_count << endl ;
 253 //    double * coord = myMesh.get_coordinates();
 254 //    os << "Coordinates :" << endl ;
 255 //    string * coordinate_name = myMesh.get_coordinates_name() ;
 256 //    string * coordinate_unit = myMesh.get_coordinates_unit() ;
 257
 258 //    for(int i=0;i<space_dimension;i++) {
 259 //      os << "  - name "<< i << " : "<< coordinate_name[i] << endl;
 260 //      os << "  - unit "<< i << " : "<< coordinate_unit[i] << endl ;
 261 //    }
 262
 263 //    for(int j=0;j<nodes_count;j++) {
 264 //      for(int i=0;i<space_dimension;i++)
 265 //        os << " coord["<< i <<","<< j << "] : "<< coord[j+i*nodes_count] <<" ," ;
 266 //      os << endl;
 267 //    }
 268
 269 //    int cells_types_count = myMesh.get_cells_types_count() ;
 270 //    os << "cells_types_count : " << cells_types_count << endl ;
 271 //    for(int i=1;i<cells_types_count;i++) {
 272 //      os << "cell type : " << myMesh.get_cells_type(i) << endl ;
 273 //      os << "  - cells count : " << myMesh.get_cells_count(i) << endl ;
 274 //      int *conectivity = myMesh.get_nodal_connectivity(i) ;
 275 //      int nodes_cells_count = myMesh.get_cells_type(i).get_nodes_count() ;
 276 //      for(int j=0;j<myMesh.get_cells_count(i);j++) {
 277 //        os << "    " ;
 278 //        for(int k=0;k<nodes_cells_count;k++)
 279 //      os <<conectivity[j*nodes_cells_count+k] << " " ;
 280 //        os << endl ;
 281 //      }
 282 //    }
 283
 284 //    int nodes_families_count = myMesh.get_nodes_families_count() ;
 285 //    os << "nodes_families_count : " << nodes_families_count << endl ;
 286 //    vector<FamilyNode*> nodes_Families = myMesh.get_nodes_Families() ;
 287 //    for(int i=0;i<nodes_families_count;i++)
 288 //      os << "  - Famile "<<i<<" : "<< (FamilyNode &) (*nodes_Families[i]) << endl ;
 289 //    os << endl ;
 290
 291 //    int cells_families_count = myMesh.get_cells_families_count() ;
 292 //    os << "cells_families_count : " << cells_families_count << endl ;
 293 //    vector<FamilyCell*> cells_Families = myMesh.get_cells_Families() ;
 294 //    for(int i=0;i<cells_families_count;i++)
 295 //      os << "  - Famile "<<i<<" : "<< (*cells_Families[i]) << endl ;
 296 //    os << endl ;
 297
 298   return os ;
 299 }
 300
 301 /*!
 302   Get global number of element which have same connectivity than connectivity argument
 303 */
 304 int MESH::getElementNumber(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type, int * connectivity)
 305 {
 306   const char* LOC="MESH::getElementNumber " ;
 307   BEGIN_OF(LOC) ;
 308
 309   int numberOfValue ; // size of connectivity array
 310   CELLMODEL myType(Type) ;
 311   if (ConnectivityType==MED_DESCENDING)
 312     numberOfValue = myType.getNumberOfConstituents(1) ; // faces (3D) or edges (2D)
 313   else
 314     numberOfValue = myType.getNumberOfNodes() ; // nodes
 315
 316   int * myReverseConnectivityValue = getReverseConnectivity(ConnectivityType) ;
 317   int * myReverseConnectivityIndex = getReverseConnectivityIndex(ConnectivityType) ;
 318
 319   // First node or face/edge
 320   int indexBegin = myReverseConnectivityIndex[connectivity[0]-1] ;
 321   int indexEnd = myReverseConnectivityIndex[connectivity[0]] ;
 322
 323   // list to put cells numbers
 324   list<int> cellsList ;
 325   list<int>::iterator itList ;
 326   for (int i=indexBegin; i<indexEnd; i++)
 327     cellsList.push_back(myReverseConnectivityValue[i-1]) ;
 328
 329   // Foreach node or face/edge in cell, we search cells which are in common.
 330   // At the end, we must have only one !
 331
 332   for (int i=1; i<numberOfValue; i++) {
 333     int connectivity_i = connectivity[i] ;
 334     for (itList=cellsList.begin();itList!=cellsList.end();itList++) {
 335       bool find = false ;
 336       for (int j=myReverseConnectivityIndex[connectivity_i-1]; j<myReverseConnectivityIndex[connectivity_i]; j++) {
 337         if ((*itList)==myReverseConnectivityValue[j-1]) {
 338           find=true ;
 339           break ;
 340         }
 341       }
 342       if (!find) {
 343         itList=cellsList.erase(itList);
 344         itList--; // well : rigth if itList = cellsList.begin() ??
 345       }
 346     }
 347   }
 348
 349   if (cellsList.size()>1) // we have more than one cell
 350     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Found more than one element !")) ;
 351
 352   if (cellsList.size()==0)
 353     return -1;
 354
 355   END_OF(LOC);
 356
 357   return cellsList.front() ;
 358 }
 359
 360
 361 /*!
 362   Return a support which reference all elements on the boundary of mesh.
 363
 364   For instance, we could get only face in 3D and edge in 2D.
 365 */
 366 SUPPORT * MESH::getBoundaryElements(medEntityMesh Entity) throw (MEDEXCEPTION)
 367 {
 368   const char * LOC = "MESH::getBoundaryElements : " ;
 369   BEGIN_OF(LOC) ;
 370   // some test :
 371   // actually we could only get face (in 3D) and edge (in 2D)
 372   if (_spaceDimension == 3)
 373     if (Entity != MED_FACE)
 374       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<Entity<<" !"));
 375   if (_spaceDimension == 2)
 376     if (Entity != MED_EDGE)
 377       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<Entity<<" !"));
 378
 379   // well, all rigth !
 380   SUPPORT * mySupport = new SUPPORT(this,"Boundary",Entity);
 381   //mySupport.setDescription("boundary of type ...");
 382   mySupport->setAll(false);
 383
 384
 385   int * myConnectivityValue = getReverseConnectivity(MED_DESCENDING) ;
 386   int * myConnectivityIndex = getReverseConnectivityIndex(MED_DESCENDING) ;
 387   int numberOf = getNumberOfElements(Entity,MED_ALL_ELEMENTS) ;
 388   list<int> myElementsList ;
 389   int size = 0 ;
 390   SCRUTE(numberOf) ;
 391   for (int i=0 ; i<numberOf; i++)
 392     if (myConnectivityValue[myConnectivityIndex[i]] == 0) {
 393       SCRUTE(i+1) ;
 394       myElementsList.push_back(i+1) ;
 395       size++ ;
 396     }
 397   SCRUTE(size) ;
 398   // Well, we must know how many geometric type we have found
 399   int * myListArray = new int[size] ;
 400   int id = 0 ;
 401   list<int>::iterator myElementsListIt ;
 402   for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
 403     myListArray[id]=(*myElementsListIt) ;
 404     SCRUTE(id);
 405     SCRUTE(myListArray[id]);
 406     id ++ ;
 407   }
 408
 409   int numberOfGeometricType ;
 410   medGeometryElement* geometricType ;
 411   int * numberOfGaussPoint ;
 412   int * geometricTypeNumber ;
 413   int * numberOfEntities ;
 414   MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY() ;
 415   int * mySkyLineArrayIndex ;
 416
 417   int numberOfType = getNumberOfTypes(Entity) ;
 418   if (numberOfType == 1) { // wonderfull : it's easy !
 419     numberOfGeometricType = 1 ;
 420     geometricType = new medGeometryElement[1] ;
 421     medGeometryElement *  allType = getTypes(Entity);
 422     geometricType[0] = allType[0] ;
 423     numberOfGaussPoint = new int[1] ;
 424     numberOfGaussPoint[0] = 1 ;
 425     geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
 426     geometricTypeNumber[0] = 0 ;
 427     numberOfEntities = new int[1] ;
 428     numberOfEntities[0] = size ;
 429     mySkyLineArrayIndex = new int[2] ;
 430     mySkyLineArrayIndex[0]=1 ;
 431     mySkyLineArrayIndex[1]=1+size ;
 432   }
 433   else {// hemmm
 434     map<medGeometryElement,int> theType ;
 435     for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
 436       medGeometryElement myType = getElementType(Entity,*myElementsListIt) ;
 437       if (theType.find(myType) != theType.end() )
 438         theType[myType]+=1 ;
 439       else
 440         theType[myType]=1 ;
 441     }
 442     numberOfGeometricType = theType.size() ;
 443     geometricType = new medGeometryElement[numberOfGeometricType] ;
 444     medGeometryElement *  allType = getTypes(Entity);
 445     numberOfGaussPoint = new int[numberOfGeometricType] ;
 446     geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
 447     numberOfEntities = new int[numberOfGeometricType] ;
 448     mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
 449     int index = 0 ;
 450     mySkyLineArrayIndex[0]=1 ;
 451     map<medGeometryElement,int>::iterator theTypeIt ;
 452     for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
 453       geometricType[index] = (*theTypeIt).first ;
 454       numberOfGaussPoint[index] = 1 ;
 455       geometricTypeNumber[index] = 0 ;
 456       numberOfEntities[index] = (*theTypeIt).second ;
 457       mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfEntities[index] ;
 458       index++ ;
 459     }
 460   }
 461   mySkyLineArray->setMEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
 462
 463   mySupport->setNumberOfGeometricType(numberOfGeometricType) ;
 464   mySupport->setGeometricType(geometricType) ;
 465   mySupport->setNumberOfGaussPoint(numberOfGaussPoint) ;
 466   mySupport->setGeometricTypeNumber(geometricTypeNumber) ;
 467   mySupport->setNumberOfEntities(numberOfEntities) ;
 468   mySupport->setTotalNumberOfEntities(size) ;
 469   mySupport->setNumber(mySkyLineArray) ;
 470
 471   END_OF(LOC) ;
 472   return mySupport ;
 473 }
 474
 475 FIELD<double>* MESH::getVolume(const SUPPORT * Support) throw (MEDEXCEPTION)
 476 {
 477   // Support must be on 3D elements
 478   MESSAGE("MESH::getVolume(SUPPORT*)");
 479
 480   // Make sure that the MESH class is the same as the MESH class attribut
 481   // in the class Support
 482   MESH* myMesh = Support->getMesh();
 483   if (this != myMesh)
 484     throw MEDEXCEPTION("MESH::getVolume(SUPPORT*) no compatibility between *this and SUPPORT::_mesh !");
 485
 486   int dim_space = getSpaceDimension();
 487   medEntityMesh support_entity = Support->getEntity();
 488   bool onAll = Support->isOnAllElements();
 489
 490   int nb_type, length_values;
 491   medGeometryElement* types;
 492   int nb_entity_type;
 493   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
 494   int* global_connectivity;
 495
 496   if (onAll)
 497     {
 498       nb_type = myMesh->getNumberOfTypes(support_entity);
 499       length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
 500       types = getTypes(support_entity);
 501     }
 502   else
 503     {
 504       nb_type = Support->getNumberOfTypes();
 505       length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
 506       types = Support->getTypes();
 507     }
 508
 509   int index;
 510   FIELD<double>* Volume;
 511
 512   Volume = new FIELD<double>::FIELD(Support,1);
 513   //  double *volume = new double [length_values];
 514   Volume->setName("VOLUME");
 515   Volume->setDescription("cells volume");
 516   Volume->setComponentName(1,"volume");
 517   Volume->setComponentDescription(1,"desc-comp");
 518
 519   /*  string MEDComponentUnit(MED_TAILLE_PNOM,' ');*/
 520
 521   string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1]+"*"+myMesh->getCoordinatesUnits()[2];
 522
 523   Volume->setMEDComponentUnit(1,MEDComponentUnit);
 524
 525   Volume->setValueType(MED_REEL64);
 526
 527   Volume->setIterationNumber(0);
 528   Volume->setOrderNumber(0);
 529   Volume->setTime(0.0);
 530
 531   double *volume = Volume->getValue(MED_FULL_INTERLACE);
 532
 533   index = 0;
 534   const double * coord = getCoordinates(MED_FULL_INTERLACE);
 535
 536   for (int i=0;i<nb_type;i++)
 537     {
 538       medGeometryElement type = types[i] ;
 539       double xvolume;
 540
 541       if (onAll)
 542         {
 543           nb_entity_type = getNumberOfElements(support_entity,type);
 544           global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
 545         }
 546       else
 547         {
 548           nb_entity_type = Support->getNumberOfElements(type);
 549
 550           int * supp_number = Support->getNumber(type);
 551           int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
 552           int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
 553           global_connectivity = new int[(type%100)*nb_entity_type];
 554
 555           for (int k_type = 0; k_type<nb_entity_type; k_type++) {
 556             for (int j_ent = 0; j_ent<(type%100); j_ent++) {
 557               global_connectivity[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
 558             }
 559           }
 560         }
 561
 562       switch (type)
 563         {
 564         case MED_TETRA4 : case MED_TETRA10 :
 565           {
 566             for (int tetra=0;tetra<nb_entity_type;tetra++)
 567               {
 568                 int tetra_index = (type%100)*tetra;
 569
 570                 int N1 = global_connectivity[tetra_index]-1;
 571                 int N2 = global_connectivity[tetra_index+1]-1;
 572                 int N3 = global_connectivity[tetra_index+2]-1;
 573                 int N4 = global_connectivity[tetra_index+3]-1;
 574
 575                 double x1 = coord[dim_space*N1];
 576                 double x2 = coord[dim_space*N2];
 577                 double x3 = coord[dim_space*N3];
 578                 double x4 = coord[dim_space*N4];
 579
 580                 double y1 = coord[dim_space*N1+1];
 581                 double y2 = coord[dim_space*N2+1];
 582                 double y3 = coord[dim_space*N3+1];
 583                 double y4 = coord[dim_space*N4+1];
 584
 585                 double z1 = coord[dim_space*N1+2];
 586                 double z2 = coord[dim_space*N2+2];
 587                 double z3 = coord[dim_space*N3+2];
 588                 double z4 = coord[dim_space*N4+2];
 589
 590                 xvolume = ((x3-x1)*((y2-y1)*(z4-z1) - (z2-z1)*(y4-y1)) -
 591                            (x2-x1)*((y3-y1)*(z4-z1) - (z3-z1)*(y4-y1)) +
 592                            (x4-x1)*((y3-y1)*(z2-z1) - (z3-z1)*(y2-y1)))/6.0;
 593
 594                 volume[index] = xvolume ;
 595                 index++;
 596               }
 597             break;
 598           }
 599         case MED_PYRA5 : case MED_PYRA13 :
 600           {
 601             for (int pyra=0;pyra<nb_entity_type;pyra++)
 602               {
 603                 int pyra_index = (type%100)*pyra;
 604
 605                 int N1 = global_connectivity[pyra_index]-1;
 606                 int N2 = global_connectivity[pyra_index+1]-1;
 607                 int N3 = global_connectivity[pyra_index+2]-1;
 608                 int N4 = global_connectivity[pyra_index+3]-1;
 609                 int N5 = global_connectivity[pyra_index+4]-1;
 610
 611                 double x1 = coord[dim_space*N1];
 612                 double x2 = coord[dim_space*N2];
 613                 double x3 = coord[dim_space*N3];
 614                 double x4 = coord[dim_space*N4];
 615                 double x5 = coord[dim_space*N5];
 616
 617                 double y1 = coord[dim_space*N1+1];
 618                 double y2 = coord[dim_space*N2+1];
 619                 double y3 = coord[dim_space*N3+1];
 620                 double y4 = coord[dim_space*N4+1];
 621                 double y5 = coord[dim_space*N5+1];
 622
 623                 double z1 = coord[dim_space*N1+2];
 624                 double z2 = coord[dim_space*N2+2];
 625                 double z3 = coord[dim_space*N3+2];
 626                 double z4 = coord[dim_space*N4+2];
 627                 double z5 = coord[dim_space*N5+2];
 628
 629                 xvolume = (((x3-x1)*((y2-y1)*(z5-z1) - (z2-z1)*(y5-y1)) -
 630                             (x2-x1)*((y3-y1)*(z5-z1) - (z3-z1)*(y5-y1)) +
 631                             (x5-x1)*((y3-y1)*(z2-z1) - (z3-z1)*(y2-y1))) +
 632                            ((x4-x1)*((y3-y1)*(z5-z1) - (z3-z1)*(y5-y1)) -
 633                             (x3-x1)*((y4-y1)*(z5-z1) - (z4-z1)*(y5-y1)) +
 634                             (x5-x1)*((y4-y1)*(z3-z1) - (z4-z1)*(y3-y1)))
 635                            )/6.0;
 636
 637                 volume[index] = xvolume ;
 638                 index = index++;
 639               }
 640             break;
 641           }
 642         case MED_PENTA6 : case MED_PENTA15 :
 643           {
 644             for (int penta=0;penta<nb_entity_type;penta++)
 645               {
 646                 int penta_index = (type%100)*penta;
 647
 648                 int N1 = global_connectivity[penta_index]-1;
 649                 int N2 = global_connectivity[penta_index+1]-1;
 650                 int N3 = global_connectivity[penta_index+2]-1;
 651                 int N4 = global_connectivity[penta_index+3]-1;
 652                 int N5 = global_connectivity[penta_index+4]-1;
 653                 int N6 = global_connectivity[penta_index+5]-1;
 654
 655                 double x1 = coord[dim_space*N1];
 656                 double x2 = coord[dim_space*N2];
 657                 double x3 = coord[dim_space*N3];
 658                 double x4 = coord[dim_space*N4];
 659                 double x5 = coord[dim_space*N5];
 660                 double x6 = coord[dim_space*N6];
 661
 662                 double y1 = coord[dim_space*N1+1];
 663                 double y2 = coord[dim_space*N2+1];
 664                 double y3 = coord[dim_space*N3+1];
 665                 double y4 = coord[dim_space*N4+1];
 666                 double y5 = coord[dim_space*N5+1];
 667                 double y6 = coord[dim_space*N6+1];
 668
 669                 double z1 = coord[dim_space*N1+2];
 670                 double z2 = coord[dim_space*N2+2];
 671                 double z3 = coord[dim_space*N3+2];
 672                 double z4 = coord[dim_space*N4+2];
 673                 double z5 = coord[dim_space*N5+2];
 674                 double z6 = coord[dim_space*N6+2];
 675
 676                 double a1 = (x2-x3)/2.0, a2 = (y2-y3)/2.0, a3 = (z2-z3)/2.0;
 677                 double b1 = (x5-x6)/2.0, b2 = (y5-y6)/2.0, b3 = (z5-z6)/2.0;
 678                 double c1 = (x4-x1)/2.0, c2 = (y4-y1)/2.0, c3 = (z4-z1)/2.0;
 679                 double d1 = (x5-x2)/2.0, d2 = (y5-y2)/2.0, d3 = (z5-z2)/2.0;
 680                 double e1 = (x6-x3)/2.0, e2 = (y6-y3)/2.0, e3 = (z6-z3)/2.0;
 681                 double f1 = (x1-x3)/2.0, f2 = (y1-y3)/2.0, f3 = (z1-z3)/2.0;
 682                 double h1 = (x4-x6)/2.0, h2 = (y4-y6)/2.0, h3 = (z4-z6)/2.0;
 683
 684                 double A = a1*c2*f3 - a1*c3*f2 - a2*c1*f3 + a2*c3*f1 +
 685                   a3*c1*f2 - a3*c2*f1;
 686                 double B = b1*c2*h3 - b1*c3*h2 - b2*c1*h3 + b2*c3*h1 +
 687                   b3*c1*h2 - b3*c2*h1;
 688                 double C = (a1*c2*h3 + b1*c2*f3) - (a1*c3*h2 + b1*c3*f2) -
 689                   (a2*c1*h3 + b2*c1*f3) + (a2*c3*h1 + b2*c3*f1) +
 690                   (a3*c1*h2 + b3*c1*f2) - (a3*c2*h1 + b3*c2*f1);
 691                 double D = a1*d2*f3 - a1*d3*f2 - a2*d1*f3 + a2*d3*f1 +
 692                   a3*d1*f2 - a3*d2*f1;
 693                 double E = b1*d2*h3 - b1*d3*h2 - b2*d1*h3 + b2*d3*h1 +
 694                   b3*d1*h2 - b3*d2*h1;
 695                 double F = (a1*d2*h3 + b1*d2*f3) - (a1*d3*h2 + b1*d3*f2) -
 696                   (a2*d1*h3 + b2*d1*f3) + (a2*d3*h1 + b2*d3*f1) +
 697                   (a3*d1*h2 + b3*d1*f2) - (a3*d2*h1 + b3*d2*f1);
 698                 double G = a1*e2*f3 - a1*e3*f2 - a2*e1*f3 + a2*e3*f1 +
 699                   a3*e1*f2 - a3*e2*f1;
 700                 double H = b1*e2*h3 - b1*e3*h2 - b2*e1*h3 + b2*e3*h1 +
 701                   b3*e1*h2 - b3*e2*h1;
 702                 double P = (a1*e2*h3 + b1*e2*f3) - (a1*e3*h2 + b1*e3*f2) -
 703                   (a2*e1*h3 + b2*e1*f3) + (a2*e3*h1 + b2*e3*f1) +
 704                   (a3*e1*h2 + b3*e1*f2) - (a3*e2*h1 + b3*e2*f1);
 705
 706                 xvolume = -2.0*(2.0*(A + B + D + E + G + H) + C + F + P)/9.0;
 707
 708                 volume[index] = xvolume ;
 709                 index++;
 710               }
 711             break;
 712           }
 713         case MED_HEXA8 : case MED_HEXA20 :
 714           {
 715             for (int hexa=0;hexa<nb_entity_type;hexa++)
 716               {
 717                 int hexa_index = (type%100)*hexa;
 718
 719                 int N1 = global_connectivity[hexa_index]-1;
 720                 int N2 = global_connectivity[hexa_index+1]-1;
 721                 int N3 = global_connectivity[hexa_index+2]-1;
 722                 int N4 = global_connectivity[hexa_index+3]-1;
 723                 int N5 = global_connectivity[hexa_index+4]-1;
 724                 int N6 = global_connectivity[hexa_index+5]-1;
 725                 int N7 = global_connectivity[hexa_index+6]-1;
 726                 int N8 = global_connectivity[hexa_index+7]-1;
 727
 728                 double x1 = coord[dim_space*N1];
 729                 double x2 = coord[dim_space*N2];
 730                 double x3 = coord[dim_space*N3];
 731                 double x4 = coord[dim_space*N4];
 732                 double x5 = coord[dim_space*N5];
 733                 double x6 = coord[dim_space*N6];
 734                 double x7 = coord[dim_space*N7];
 735                 double x8 = coord[dim_space*N8];
 736
 737                 double y1 = coord[dim_space*N1+1];
 738                 double y2 = coord[dim_space*N2+1];
 739                 double y3 = coord[dim_space*N3+1];
 740                 double y4 = coord[dim_space*N4+1];
 741                 double y5 = coord[dim_space*N5+1];
 742                 double y6 = coord[dim_space*N6+1];
 743                 double y7 = coord[dim_space*N7+1];
 744                 double y8 = coord[dim_space*N8+1];
 745
 746                 double z1 = coord[dim_space*N1+2];
 747                 double z2 = coord[dim_space*N2+2];
 748                 double z3 = coord[dim_space*N3+2];
 749                 double z4 = coord[dim_space*N4+2];
 750                 double z5 = coord[dim_space*N5+2];
 751                 double z6 = coord[dim_space*N6+2];
 752                 double z7 = coord[dim_space*N7+2];
 753                 double z8 = coord[dim_space*N8+2];
 754
 755                 double a1 = (x3-x4)/8.0, a2 = (y3-y4)/8.0, a3 = (z3-z4)/8.0;
 756                 double b1 = (x2-x1)/8.0, b2 = (y2-y1)/8.0, b3 = (z2-z1)/8.0;
 757                 double c1 = (x7-x8)/8.0, c2 = (y7-y8)/8.0, c3 = (z7-z8)/8.0;
 758                 double d1 = (x6-x5)/8.0, d2 = (y6-y5)/8.0, d3 = (z6-z5)/8.0;
 759                 double e1 = (x3-x2)/8.0, e2 = (y3-y2)/8.0, e3 = (z3-z2)/8.0;
 760                 double f1 = (x4-x1)/8.0, f2 = (y4-y1)/8.0, f3 = (z4-z1)/8.0;
 761                 double h1 = (x7-x6)/8.0, h2 = (y7-y6)/8.0, h3 = (z7-z6)/8.0;
 762                 double p1 = (x8-x5)/8.0, p2 = (y8-y5)/8.0, p3 = (z8-z5)/8.0;
 763                 double q1 = (x3-x7)/8.0, q2 = (y3-y7)/8.0, q3 = (z3-z7)/8.0;
 764                 double r1 = (x4-x8)/8.0, r2 = (y4-y8)/8.0, r3 = (z4-z8)/8.0;
 765                 double s1 = (x2-x6)/8.0, s2 = (y2-y6)/8.0, s3 = (z2-z6)/8.0;
 766                 double t1 = (x1-x5)/8.0, t2 = (y1-y5)/8.0, t3 = (z1-z5)/8.0;
 767
 768                 double A = a1*e2*q3 - a1*e3*q2 - a2*e1*q3 + a2*e3*q1 +
 769                   a3*e1*q2 - a3*e2*q1;
 770                 double B = c1*h2*q3 - c1*h3*q2 - c2*h1*q3 + c2*h3*q1 +
 771                   c3*h1*q2 - c3*h2*q1;
 772                 double C = (a1*h2 + c1*e2)*q3 - (a1*h3 + c1*e3)*q2 -
 773                   (a2*h1 + c2*e1)*q3 + (a2*h3 + c2*e3)*q1 +
 774                   (a3*h1 + c3*e1)*q2 - (a3*h2 + c3*e2)*q1;
 775                 double D = b1*e2*s3 - b1*e3*s2 - b2*e1*s3 + b2*e3*s1 +
 776                   b3*e1*s2 - b3*e2*s1;
 777                 double E = d1*h2*s3 - d1*h3*s2 - d2*h1*s3 + d2*h3*s1 +
 778                   d3*h1*s2 - d3*h2*s1;
 779                 double F = (b1*h2 + d1*e2)*s3 - (b1*h3 + d1*e3)*s2 -
 780                   (b2*h1 + d2*e1)*s3 + (b2*h3 + d2*e3)*s1 +
 781                   (b3*h1 + d3*e1)*s2 - (b3*h2 + d3*e2)*s1;
 782                 double G = (a1*e2*s3 + b1*e2*q3) - (a1*e3*s2 + b1*e3*q2) -
 783                   (a2*e1*s3 + b2*e1*q3) + (a2*e3*s1 + b2*e3*q1) +
 784                   (a3*e1*s2 + b3*e1*q2) - (a3*e2*s1 + b3*e2*q1);
 785                 double H = (c1*h2*s3 + d1*h2*q3) - (c1*h3*s2 + d1*h3*q2) -
 786                   (c2*h1*s3 + d2*h1*q3) + (c2*h3*s1 + d2*h3*q1) +
 787                   (c3*h1*s2 + d3*h1*q2) - (c3*h2*s1 + d3*h2*q1);
 788                 double I = ((a1*h2 + c1*e2)*s3 + (b1*h2 + d1*e2)*q3) -
 789                   ((a1*h3 + c1*e3)*s2 + (b1*h3 + d1*e3)*q2) -
 790                   ((a2*h1 + c2*e1)*s3 + (b2*h1 + d2*e1)*q3) +
 791                   ((a2*h3 + c2*e3)*s1 + (b2*h3 + d2*e3)*q1) +
 792                   ((a3*h1 + c3*e1)*s2 + (b3*h1 + d3*e1)*q2) -
 793                   ((a3*h2 + c3*e2)*s1 + (b3*h2 + d3*e2)*q1);
 794                 double J = a1*f2*r3 - a1*f3*r2 - a2*f1*r3 + a2*f3*r1 +
 795                   a3*f1*r2 - a3*f2*r1;
 796                 double K = c1*p2*r3 - c1*p3*r2 - c2*p1*r3 + c2*p3*r1 +
 797                   c3*p1*r2 - c3*p2*r1;
 798                 double L = (a1*p2 + c1*f2)*r3 - (a1*p3 + c1*f3)*r2 -
 799                   (a2*p1 + c2*f1)*r3 + (a2*p3 + c2*f3)*r1 +
 800                   (a3*p1 + c3*f1)*r2 - (a3*p2 + c3*f2)*r1;
 801                 double M = b1*f2*t3 - b1*f3*t2 - b2*f1*t3 + b2*f3*t1 +
 802                   b3*f1*t2 - b3*f2*t1;
 803                 double N = d1*p2*t3 - d1*p3*t2 - d2*p1*t3 + d2*p3*t1 +
 804                   d3*p1*t2 - d3*p2*t1;
 805                 double O = (b1*p2 + d1*f2)*t3 - (b1*p3 + d1*f3)*t2 -
 806                   (b2*p1 + d2*f1)*t3 + (b2*p3 + d2*f3)*t1 +
 807                   (b3*p1 + d3*f1)*t2 - (b3*p2 + d3*f2)*t1;
 808                 double P = (a1*f2*t3 + b1*f2*r3) - (a1*f3*t2 + b1*f3*r2) -
 809                   (a2*f1*t3 + b2*f1*r3) + (a2*f3*t1 + b2*f3*r1) +
 810                   (a3*f1*t2 + b3*f1*r2) - (a3*f2*t1 + b3*f2*r1);
 811                 double Q = (c1*p2*t3 + d1*p2*r3) - (c1*p3*t2 + d1*p3*r2) -
 812                   (c2*p1*t3 + d2*p1*r3) + (c2*p3*t1 + d2*p3*r1) +
 813                   (c3*p1*t2 + d3*p1*r2) - (c3*p2*t1 + d3*p2*r1);
 814                 double R = ((a1*p2 + c1*f2)*t3 + (b1*p2 + d1*f2)*r3) -
 815                   ((a1*p3 + c1*f3)*t2 + (b1*p3 + d1*f3)*r2) -
 816                   ((a2*p1 + c2*f1)*t3 + (b2*p1 + d2*f1)*r3) +
 817                   ((a2*p3 + c2*f3)*t1 + (b2*p3 + d2*f3)*r1) +
 818                   ((a3*p1 + c3*f1)*t2 + (b3*p1 + d3*f1)*r2) -
 819                   ((a3*p2 + c3*f2)*t1 + (b3*p2 + d3*f2)*r1);
 820                 double S = (a1*e2*r3 + a1*f2*q3) - (a1*e3*r2 + a1*f3*q2) -
 821                   (a2*e1*r3 + a2*f1*q3) + (a2*e3*r1 + a2*f3*q1) +
 822                   (a3*e1*r2 + a3*f1*q2) - (a3*e2*r1 + a3*f2*q1);
 823                 double T = (c1*h2*r3 + c1*p2*q3) - (c1*h3*r2 + c1*p3*q2) -
 824                   (c2*h1*r3 + c2*p1*q3) + (c2*h3*r1 + c2*p3*q1) +
 825                   (c3*h1*r2 + c3*p1*q2) - (c3*h2*r1 + c3*p2*q1);
 826                 double U = ((a1*h2 + c1*e2)*r3 + (a1*p2 + c1*f2)*q3) -
 827                   ((a1*h3 + c1*e3)*r2 + (a1*p3 + c1*f3)*q2) -
 828                   ((a2*h1 + c2*e1)*r3 + (a2*p1 + c2*f1)*q3) +
 829                   ((a2*h3 + c2*e3)*r1 + (a2*p3 + c2*f3)*q1) +
 830                   ((a3*h1 + c3*e1)*r2 + (a3*p1 + c3*f1)*q2) -
 831                   ((a3*h2 + c3*e2)*r1 + (a3*p2 + c3*f2)*q1);
 832                 double V = (b1*e2*t3 + b1*f2*s3) - (b1*e3*t2 + b1*f3*s2) -
 833                   (b2*e1*t3 + b2*f1*s3) + (b2*e3*t1 + b2*f3*s1) +
 834                   (b3*e1*t2 + b3*f1*s2) - (b3*e2*t1 + b3*f2*s1);
 835                 double W = (d1*h2*t3 + d1*p2*s3) - (d1*h3*t2 + d1*p3*s2) -
 836                   (d2*h1*t3 + d2*p1*s3) + (d2*h3*t1 + d2*p3*s1) +
 837                   (d3*h1*t2 + d3*p1*s2) - (d3*h2*t1 + d3*p2*s1);
 838                 double X = ((b1*h2 + d1*e2)*t3 + (b1*p2 + d1*f2)*s3) -
 839                   ((b1*h3 + d1*e3)*t2 + (b1*p3 + d1*f3)*s2) -
 840                   ((b2*h1 + d2*e1)*t3 + (b2*p1 + d2*f1)*s3) +
 841                   ((b2*h3 + d2*e3)*t1 + (b2*p3 + d2*f3)*s1) +
 842                   ((b3*h1 + d3*e1)*t2 + (b3*p1 + d3*f1)*s2) -
 843                   ((b3*h2 + d3*e2)*t1 + (b3*p2 + d3*f2)*s1);
 844                 double Y = (a1*e2*t3 + a1*f2*s3 + b1*e2*r3 + b1*f2*q3) -
 845                   (a1*e3*t2 + a1*f3*s2 + b1*e3*r2 + b1*f3*q2) -
 846                   (a2*e1*t3 + a2*f1*s3 + b2*e1*r3 + b2*f1*q3) +
 847                   (a2*e3*t1 + a2*f3*s1 + b2*e3*r1 + b2*f3*q1) +
 848                   (a3*e1*t2 + a3*f1*s2 + b3*e1*r2 + b3*f1*q2) -
 849                   (a3*e2*t1 + a3*f2*s1 + b3*e2*r1 + b3*f2*q1);
 850                 double Z = (c1*h2*t3 + c1*p2*s3 + d1*h2*r3 + d1*p2*q3) -
 851                   (c1*h3*t2 + c1*p3*s2 + d1*h3*r2 + d1*p3*q2) -
 852                   (c2*h1*t3 + c2*p1*s3 + d2*h1*r3 + d2*p1*q3) +
 853                   (c2*h3*t1 + c2*p3*s1 + d2*h3*r1 + d2*p3*q1) +
 854                   (c3*h1*t2 + c3*p1*s2 + d3*h1*r2 + d3*p1*q2) -
 855                   (c3*h2*t1 + c3*p2*s1 + d3*h2*r1 + d3*p2*q1);
 856                 double AA = ((a1*h2 + c1*e2)*t3 + (a1*p2 + c1*f2)*s3 +
 857                              (b1*h2 + d1*e2)*r3 + (b1*p2 + d1*f2)*q3) -
 858                   ((a1*h3 + c1*e3)*t2 + (a1*p3 + c1*f3)*s2 +
 859                    (b1*h3 + d1*e3)*r2 + (b1*p3 + d1*f3)*q2) -
 860                   ((a2*h1 + c2*e1)*t3 + (a2*p1 + c2*f1)*s3 +
 861                    (b2*h1 + d2*e1)*r3 + (b2*p1 + d2*f1)*q3) +
 862                   ((a2*h3 + c2*e3)*t1 + (a2*p3 + c2*f3)*s1 +
 863                    (b2*h3 + d2*e3)*r1 + (b2*p3 + d2*f3)*q1) +
 864                   ((a3*h1 + c3*e1)*t2 + (a3*p1 + c3*f1)*s2 +
 865                    (b3*h1 + d3*e1)*r2 + (b3*p1 + d3*f1)*q2) -
 866                   ((a3*h2 + c3*e2)*t1 + (a3*p2 + c3*f2)*s1 +
 867                    (b3*h2 + d3*e2)*r1 + (b3*p2 + d3*f2)*q1);
 868
 869                 xvolume = 64.0*(8.0*(A + B + D + E + J + K + M + N) +
 870                                 4.0*(C + F + G + H + L + O + P + Q + S + T +
 871                                      V + W) + 2.0*(I + R + U + X + Y + Z) +
 872                                 AA)/27.0;
 873
 874                 volume[index] = xvolume ;
 875                 index++;
 876               }
 877             break;
 878           }
 879         default :
 880           throw MEDEXCEPTION("MESH::getVolume(SUPPORT*) Bad Support to get Volumes on it !");
 881           break;
 882         }
 883     }
 884
 885   return Volume;
 886 }
 887
 888 FIELD<double>* MESH::getArea(const SUPPORT * Support) throw (MEDEXCEPTION)
 889 {
 890   // Support must be on 2D elements
 891   MESSAGE("MESH::getArea(SUPPORT*)");
 892
 893   // Make sure that the MESH class is the same as the MESH class attribut
 894   // in the class Support
 895   MESH* myMesh = Support->getMesh();
 896   if (this != myMesh)
 897     throw MEDEXCEPTION("MESH::getArea(SUPPORT*) no compatibility between *this and SUPPORT::_mesh !");
 898
 899   int dim_space = getSpaceDimension();
 900   medEntityMesh support_entity = Support->getEntity();
 901   bool onAll = Support->isOnAllElements();
 902
 903   int nb_type, length_values;
 904   medGeometryElement* types;
 905   int nb_entity_type;
 906   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
 907   int* global_connectivity;
 908
 909   if (onAll)
 910     {
 911       nb_type = myMesh->getNumberOfTypes(support_entity);
 912       length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
 913       types = getTypes(support_entity);
 914     }
 915   else
 916     {
 917       nb_type = Support->getNumberOfTypes();
 918       length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
 919       types = Support->getTypes();
 920     }
 921
 922   int index;
 923   FIELD<double>* Area;
 924
 925   Area = new FIELD<double>::FIELD(Support,1);
 926   Area->setName("AREA");
 927   Area->setDescription("cells or faces area");
 928
 929   Area->setComponentName(1,"area");
 930   Area->setComponentDescription(1,"desc-comp");
 931
 932   /*  string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1];*/
 933
 934   string MEDComponentUnit="trucmuch";
 935
 936   Area->setMEDComponentUnit(1,MEDComponentUnit);
 937
 938   Area->setValueType(MED_REEL64);
 939
 940   Area->setIterationNumber(0);
 941   Area->setOrderNumber(0);
 942   Area->setTime(0.0);
 943
 944   double *area = new double[length_values];
 945   //double *area = Area->getValue(MED_FULL_INTERLACE);
 946
 947   const double * coord = getCoordinates(MED_FULL_INTERLACE);
 948   index = 0;
 949
 950   for (int i=0;i<nb_type;i++)
 951     {
 952       medGeometryElement type = types[i] ;
 953       double xarea;
 954
 955       if (onAll)
 956         {
 957           nb_entity_type = getNumberOfElements(support_entity,type);
 958           global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
 959         }
 960       else
 961         {
 962           nb_entity_type = Support->getNumberOfElements(type);
 963
 964           int * supp_number = Support->getNumber(type);
 965           int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
 966           int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
 967           global_connectivity = new int[(type%100)*nb_entity_type];
 968
 969           for (int k_type = 0; k_type<nb_entity_type; k_type++) {
 970             for (int j_ent = 0; j_ent<(type%100); j_ent++) {
 971               global_connectivity[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
 972             }
 973           }
 974         }
 975
 976       switch (type)
 977         {
 978         case MED_TRIA3 : case MED_TRIA6 :
 979           {
 980             for (int tria=0;tria<nb_entity_type;tria++)
 981               {
 982                 int tria_index = (type%100)*tria;
 983
 984                 int N1 = global_connectivity[tria_index]-1;
 985                 int N2 = global_connectivity[tria_index+1]-1;
 986                 int N3 = global_connectivity[tria_index+2]-1;
 987
 988                 double x1 = coord[dim_space*N1];
 989                 double x2 = coord[dim_space*N2];
 990                 double x3 = coord[dim_space*N3];
 991
 992                 double y1 = coord[dim_space*N1+1];
 993                 double y2 = coord[dim_space*N2+1];
 994                 double y3 = coord[dim_space*N3+1];
 995
 996                 if (dim_space==2)
 997                   {
 998                     xarea = - ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))/2.0;
 999                   }
1000                 else
1001                   {
1002                     double z1 = coord[dim_space*N1+2];
1003                     double z2 = coord[dim_space*N2+2];
1004                     double z3 = coord[dim_space*N3+2];
1005
1006                     xarea = sqrt(((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1))*
1007                                  ((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1)) +
1008                                  ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1))*
1009                                  ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1)) +
1010                                  ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))*
1011                                  ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1)))/2.0;
1012                   }
1013
1014                 area[index] = xarea ;
1015                 index++;
1016               }
1017             break;
1018           }
1019         case MED_QUAD4 : case MED_QUAD8 :
1020           {
1021             for (int quad=0;quad<nb_entity_type;quad++)
1022               {
1023                 int quad_index = (type%100)*quad;
1024
1025                 int N1 = global_connectivity[quad_index]-1;
1026                 int N2 = global_connectivity[quad_index+1]-1;
1027                 int N3 = global_connectivity[quad_index+2]-1;
1028                 int N4 = global_connectivity[quad_index+3]-1;
1029
1030                 double x1 = coord[dim_space*N1];
1031                 double x2 = coord[dim_space*N2];
1032                 double x3 = coord[dim_space*N3];
1033                 double x4 = coord[dim_space*N4];
1034
1035                 double y1 = coord[dim_space*N1+1];
1036                 double y2 = coord[dim_space*N2+1];
1037                 double y3 = coord[dim_space*N3+1];
1038                 double y4 = coord[dim_space*N4+1];
1039
1040                 if (dim_space==2)
1041                   {
1042                     double a1 = (x2-x1)/4.0, a2 = (y2-y1)/4.0;
1043                     double b1 = (x3-x4)/4.0, b2 = (y3-y4)/4.0;
1044                     double c1 = (x3-x2)/4.0, c2 = (y3-y2)/4.0;
1045                     double d1 = (x4-x1)/4.0, d2 = (y4-y1)/4.0;
1046
1047                     xarea = - 4.0*(b1*c2 - c1*b2 + a1*c2 - c1*a2 + b1*d2 -
1048                                  d1*b2 + a1*d2 - d1*a2);
1049                   }
1050                 else
1051                   {
1052                     double z1 = coord[dim_space*N1+2];
1053                     double z2 = coord[dim_space*N2+2];
1054                     double z3 = coord[dim_space*N3+2];
1055                     double z4 = coord[dim_space*N4+2];
1056
1057                     xarea = (sqrt(((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1))*
1058                                   ((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1)) +
1059                                   ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1))*
1060                                   ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1)) +
1061                                   ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))*
1062                                   ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))) +
1063                              sqrt(((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3))*
1064                                   ((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3)) +
1065                                   ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3))*
1066                                   ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3)) +
1067                                   ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))*
1068                                   ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))))/2.0;
1069                   }
1070
1071                 area[index] = xarea ;
1072                 index++;
1073               }
1074             break;
1075           }
1076         default :
1077           throw MEDEXCEPTION("MESH::getArea(SUPPORT*) Bad Support to get Areas on it !");
1078           break;
1079         }
1080     }
1081
1082   Area->setValue(MED_FULL_INTERLACE,area);
1083
1084   return Area;
1085 }
1086
1087 FIELD<double>* MESH::getLength(const SUPPORT * Support) throw (MEDEXCEPTION)
1088 {
1089   // Support must be on 1D elements
1090   MESSAGE("MESH::getLength(SUPPORT*)");
1091
1092   // Make sure that the MESH class is the same as the MESH class attribut
1093   // in the class Support
1094   MESH* myMesh = Support->getMesh();
1095   if (this != myMesh)
1096     throw MEDEXCEPTION("MESH::getLength(SUPPORT*) no compatibility between *this and SUPPORT::_mesh !");
1097
1098   int dim_space = getSpaceDimension();
1099   medEntityMesh support_entity = Support->getEntity();
1100   bool onAll = Support->isOnAllElements();
1101
1102   int nb_type, length_values;
1103   medGeometryElement* types;
1104   int nb_entity_type;
1105   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1106   int* global_connectivity;
1107
1108   if (onAll)
1109     {
1110       nb_type = myMesh->getNumberOfTypes(support_entity);
1111       length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1112       types = getTypes(support_entity);
1113     }
1114   else
1115     {
1116       nb_type = Support->getNumberOfTypes();
1117       length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1118       types = Support->getTypes();
1119     }
1120
1121   int index;
1122   FIELD<double>* Length;
1123
1124   Length = new FIELD<double>::FIELD(Support,1);
1125   //  double *length = new double [length_values];
1126   Length->setValueType(MED_REEL64);
1127
1128   double *length = Length->getValue(MED_FULL_INTERLACE);
1129
1130   const double * coord = getCoordinates(MED_FULL_INTERLACE);
1131   index = 0;
1132
1133   for (int i=0;i<nb_type;i++)
1134     {
1135       medGeometryElement type = types[i] ;
1136       double xlength;
1137
1138       if (onAll)
1139         {
1140           nb_entity_type = getNumberOfElements(support_entity,type);
1141           global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1142         }
1143       else
1144         {
1145           nb_entity_type = Support->getNumberOfElements(type);
1146
1147           int * supp_number = Support->getNumber(type);
1148           int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1149           int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1150           global_connectivity = new int[(type%100)*nb_entity_type];
1151
1152           for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1153             for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1154               global_connectivity[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1155             }
1156           }
1157         }
1158
1159       switch (type)
1160         {
1161         case MED_SEG2 : case MED_SEG3 :
1162           {
1163             for (int edge=0;edge<nb_entity_type;edge++)
1164               {
1165                 int edge_index = (type%100)*edge;
1166
1167                 int N1 = global_connectivity[edge_index]-1;
1168                 int N2 = global_connectivity[edge_index+1]-1;
1169
1170                 double x1 = coord[dim_space*N1];
1171                 double x2 = coord[dim_space*N2];
1172
1173                 double y1 = coord[dim_space*N1+1];
1174                 double y2 = coord[dim_space*N2+1];
1175
1176                 double z1, z2 ; z1 = 0.0 ; z2 = 0.0 ; if (dim_space==3) {
1177                   z1 = coord[dim_space*N1+2]; z2 = coord[dim_space*N2+2];}
1178
1179                 xlength =  sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) +
1180                                 (z1 - z2)*(z1 - z2));
1181
1182                 length[index] = xlength ;
1183                 index++;
1184               }
1185             break;
1186           }
1187         default :
1188           throw MEDEXCEPTION("MESH::getLength(SUPPORT*) Bad Support to get Lengths on it !");
1189           break;
1190         }
1191     }
1192
1193   return Length;
1194 }
1195
1196 FIELD<double>* MESH::getNormal(const SUPPORT * Support) throw (MEDEXCEPTION)
1197 {
1198   // Support must be on 2D or 1D elements
1199   MESSAGE("MESH::getNormal(SUPPORT*)");
1200
1201   // Make sure that the MESH class is the same as the MESH class attribut
1202   // in the class Support
1203   MESH* myMesh = Support->getMesh();
1204   if (this != myMesh)
1205     throw MEDEXCEPTION("MESH::getNormal(SUPPORT*) no compatibility between *this and SUPPORT::_mesh : pointeur problem !");
1206
1207   int dim_space = getSpaceDimension();
1208   medEntityMesh support_entity = Support->getEntity();
1209   bool onAll = Support->isOnAllElements();
1210
1211   int nb_type, length_values;
1212   medGeometryElement* types;
1213   int nb_entity_type;
1214   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1215   int* global_connectivity;
1216
1217   if (onAll)
1218     {
1219       nb_type = myMesh->getNumberOfTypes(support_entity);
1220       length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1221       types = getTypes(support_entity);
1222     }
1223   else
1224     {
1225       nb_type = Support->getNumberOfTypes();
1226       length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1227       types = Support->getTypes();
1228     }
1229
1230   int index;
1231
1232   FIELD<double>* Normal = new FIELD<double>::FIELD(Support,dim_space);
1233   Normal->setName("NORMAL");
1234   Normal->setDescription("cells or faces normal");
1235   for (int k=0;k<dim_space;k++) {
1236     int kp1 = k + 1;
1237     Normal->setComponentName(kp1,"normal");
1238     Normal->setComponentDescription(kp1,"desc-comp");
1239     Normal->setMEDComponentUnit(kp1,"unit");
1240   }
1241
1242   Normal->setValueType(MED_REEL64);
1243
1244   Normal->setIterationNumber(0);
1245   Normal->setOrderNumber(0);
1246   Normal->setTime(0.0);
1247
1248   double * normal = new double [dim_space*length_values];
1249   //double *normal = Normal->getValue(MED_FULL_INTERLACE);
1250
1251   const double * coord = getCoordinates(MED_FULL_INTERLACE);
1252   index = 0;
1253
1254   for (int i=0;i<nb_type;i++)
1255     {
1256       medGeometryElement type = types[i] ;
1257
1258       // Make sure we trying to get Normals on
1259       // TRIA3 or TRIA6 or QUAD4 QUAD8 (2D) cells on a 3D mesh
1260       // or on SEG2 or SEG3 (1D) cells on a 2D mesh
1261
1262       if ( (((type==MED_TRIA3) || (type==MED_TRIA6) ||
1263              (type==MED_QUAD4) || (type==MED_QUAD8)) &&
1264             (dim_space != 3)) || (((type==MED_SEG2) || (type==MED_SEG3)) &&
1265                                   (dim_space != 2)) )
1266         throw MEDEXCEPTION("MESH::getNormal(SUPPORT*) no compatibility between *this and SUPPORT::_mesh : dimension problem !");
1267
1268       double xnormal1, xnormal2, xnormal3 ;
1269
1270       if (onAll)
1271         {
1272           nb_entity_type = getNumberOfElements(support_entity,type);
1273           global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1274         }
1275       else
1276         {
1277           nb_entity_type = Support->getNumberOfElements(type);
1278
1279           int * supp_number = Support->getNumber(type);
1280           int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1281           int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1282           global_connectivity = new int[(type%100)*nb_entity_type];
1283
1284           for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1285             for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1286               global_connectivity[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1287             }
1288           }
1289         }
1290
1291       switch (type)
1292         {
1293         case MED_TRIA3 : case MED_TRIA6 :
1294           {
1295             for (int tria=0;tria<nb_entity_type;tria++)
1296               {
1297                 int tria_index = (type%100)*tria;
1298
1299                 int N1 = global_connectivity[tria_index]-1;
1300                 int N2 = global_connectivity[tria_index+1]-1;
1301                 int N3 = global_connectivity[tria_index+2]-1;
1302
1303                 double xarea;
1304                 double x1 = coord[dim_space*N1];
1305                 double x2 = coord[dim_space*N2];
1306                 double x3 = coord[dim_space*N3];
1307
1308                 double y1 = coord[dim_space*N1+1];
1309                 double y2 = coord[dim_space*N2+1];
1310                 double y3 = coord[dim_space*N3+1];
1311
1312                 double z1 = coord[dim_space*N1+2];
1313                 double z2 = coord[dim_space*N2+2];
1314                 double z3 = coord[dim_space*N3+2];
1315
1316                 xnormal1 = ((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1))/2.0;
1317                 xnormal2 = ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1))/2.0;
1318                 xnormal3 = ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))/2.0;
1319
1320                 normal[3*index] = xnormal1 ;
1321                 normal[3*index+1] = xnormal2 ;
1322                 normal[3*index+2] = xnormal3 ;
1323
1324                 index++;
1325               }
1326             break;
1327           }
1328         case MED_QUAD4 : case MED_QUAD8 :
1329           {
1330             for (int quad=0;quad<nb_entity_type;quad++)
1331               {
1332                 int quad_index = (type%100)*quad;
1333
1334                 int N1 = global_connectivity[quad_index]-1;
1335                 int N2 = global_connectivity[quad_index+1]-1;
1336                 int N3 = global_connectivity[quad_index+2]-1;
1337                 int N4 = global_connectivity[quad_index+3]-1;
1338
1339                 double xarea;
1340                 double x1 = coord[dim_space*N1];
1341                 double x2 = coord[dim_space*N2];
1342                 double x3 = coord[dim_space*N3];
1343                 double x4 = coord[dim_space*N4];
1344
1345                 double y1 = coord[dim_space*N1+1];
1346                 double y2 = coord[dim_space*N2+1];
1347                 double y3 = coord[dim_space*N3+1];
1348                 double y4 = coord[dim_space*N4+1];
1349
1350                 double z1 = coord[dim_space*N1+2];
1351                 double z2 = coord[dim_space*N2+2];
1352                 double z3 = coord[dim_space*N3+2];
1353                 double z4 = coord[dim_space*N4+2];
1354
1355                 xnormal1 = (y2-y1)*(z4-z1) - (y4-y1)*(z2-z1);
1356                 xnormal2 = (x4-x1)*(z2-z1) - (x2-x1)*(z4-z1);
1357                 xnormal3 = (x2-x1)*(y4-y1) - (x4-x1)*(y2-y1);
1358                 xarea = sqrt(xnormal1*xnormal1 + xnormal2*xnormal2 +
1359                              xnormal3*xnormal3);
1360
1361                 xnormal1 = xnormal1/xarea;
1362                 xnormal2 = xnormal2/xarea;
1363                 xnormal3 = xnormal3/xarea;
1364
1365                 xarea = (sqrt(((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1))*
1366                               ((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1)) +
1367                               ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1))*
1368                               ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1)) +
1369                               ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))*
1370                               ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))) +
1371                          sqrt(((y3-y2)*(z4-z3) - (y4-y3)*(z3-z2))*
1372                               ((y3-y2)*(z4-z3) - (y4-y3)*(z3-z2)) +
1373                               ((x4-x3)*(z2-z2) - (x3-x2)*(z4-z3))*
1374                               ((x4-x3)*(z4-z2) - (x3-x2)*(z4-z3)) +
1375                               ((x3-x2)*(y4-y3) - (x4-x3)*(y3-y2))*
1376                               ((x3-x2)*(y4-y3) - (x4-x3)*(y3-y2))))/2.0;
1377
1378                 xnormal1 = xnormal1*xarea;
1379                 xnormal2 = xnormal2*xarea;
1380                 xnormal3 = xnormal3*xarea;
1381
1382                 normal[3*index] = xnormal1 ;
1383                 normal[3*index+1] = xnormal2 ;
1384                 normal[3*index+2] = xnormal3 ;
1385
1386                 index++;
1387               }
1388             break;
1389           }
1390         case MED_SEG2 : case MED_SEG3 :
1391           {
1392             for (int edge=0;edge<nb_entity_type;edge++)
1393               {
1394                 int edge_index = (type%100)*edge;
1395
1396                 int N1 = global_connectivity[edge_index]-1;
1397                 int N2 = global_connectivity[edge_index+1]-1;
1398
1399                 double x1 = coord[dim_space*N1];
1400                 double x2 = coord[dim_space*N2];
1401
1402                 double y1 = coord[dim_space*N1+1];
1403                 double y2 = coord[dim_space*N2+1];
1404
1405                 xnormal1 = -(y2-y1);
1406                 xnormal2 = x2-x1;
1407
1408                 normal[2*index] = xnormal1 ;
1409                 normal[2*index+1] = xnormal2 ;
1410
1411                 index++;
1412               }
1413             break;
1414           }
1415         default :
1416           throw MEDEXCEPTION("MESH::getNormal(SUPPORT*) Bad Support to get Normals on it !");
1417           break;
1418         }
1419     }
1420
1421   Normal->setValue(MED_FULL_INTERLACE,normal);
1422
1423   return Normal;
1424 }
1425
1426 FIELD<double>* MESH::getBarycenter(const SUPPORT * Support) throw (MEDEXCEPTION)
1427 {
1428   MESSAGE("MESH::getBarycenter(SUPPORT*)");
1429
1430   // Make sure that the MESH class is the same as the MESH class attribut
1431   // in the class Support
1432   MESH* myMesh = Support->getMesh();
1433   if (this != myMesh)
1434     throw MEDEXCEPTION("MESH::getBarycenter(SUPPORT*) no compatibility between *this and SUPPORT::_mesh !");
1435
1436   int dim_space = getSpaceDimension();
1437   medEntityMesh support_entity = Support->getEntity();
1438   bool onAll = Support->isOnAllElements();
1439
1440   int nb_type, length_values;
1441   medGeometryElement* types;
1442   int nb_entity_type;
1443   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1444   int* global_connectivity;
1445
1446   if (onAll)
1447     {
1448       nb_type = myMesh->getNumberOfTypes(support_entity);
1449       length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1450       types = getTypes(support_entity);
1451     }
1452   else
1453     {
1454       nb_type = Support->getNumberOfTypes();
1455       length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1456       types = Support->getTypes();
1457     }
1458
1459   int index;
1460   FIELD<double>* Barycenter;
1461
1462   Barycenter = new FIELD<double>::FIELD(Support,dim_space);
1463   Barycenter->setName("BARYCENTER");
1464   Barycenter->setDescription("cells or faces barycenter");
1465
1466   for (int k=0;k<dim_space;k++) {
1467     int kp1 = k + 1;
1468     Barycenter->setComponentName(kp1,myMesh->getCoordinatesNames()[k]);
1469     Barycenter->setComponentDescription(kp1,"desc-comp");
1470     Barycenter->setMEDComponentUnit(kp1,myMesh->getCoordinatesUnits()[k]);
1471   }
1472
1473   Barycenter->setValueType(MED_REEL64);
1474
1475   Barycenter->setIterationNumber(0);
1476   Barycenter->setOrderNumber(0);
1477   Barycenter->setTime(0.0);
1478
1479   double *barycenter = new double [dim_space*length_values];
1480   //  double *barycenter = Barycenter->getValue(MED_FULL_INTERLACE);
1481
1482   const double * coord = getCoordinates(MED_FULL_INTERLACE);
1483   index = 0;
1484
1485   for (int i=0;i<nb_type;i++)
1486     {
1487       medGeometryElement type = types[i] ;
1488       double xbarycenter1, xbarycenter2, xbarycenter3;
1489
1490       if (onAll)
1491         {
1492           nb_entity_type = getNumberOfElements(support_entity,type);
1493           global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1494         }
1495       else
1496         {
1497           nb_entity_type = Support->getNumberOfElements(type);
1498
1499           int * supp_number = Support->getNumber(type);
1500           int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1501           int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1502           global_connectivity = new int[(type%100)*nb_entity_type];
1503
1504           for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1505             for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1506               global_connectivity[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1507             }
1508           }
1509         }
1510
1511       switch (type)
1512         {
1513         case MED_TETRA4 : case MED_TETRA10 :
1514           {
1515             for (int tetra =0;tetra<nb_entity_type;tetra++)
1516               {
1517                 int tetra_index = (type%100)*tetra;
1518
1519                 int N1 = global_connectivity[tetra_index]-1;
1520                 int N2 = global_connectivity[tetra_index+1]-1;
1521                 int N3 = global_connectivity[tetra_index+2]-1;
1522                 int N4 = global_connectivity[tetra_index+3]-1;
1523
1524                 double x1 = coord[dim_space*N1];
1525                 double x2 = coord[dim_space*N2];
1526                 double x3 = coord[dim_space*N3];
1527                 double x4 = coord[dim_space*N4];
1528
1529                 double y1 = coord[dim_space*N1+1];
1530                 double y2 = coord[dim_space*N2+1];
1531                 double y3 = coord[dim_space*N3+1];
1532                 double y4 = coord[dim_space*N4+1];
1533
1534                 double z1 = coord[dim_space*N1+2];
1535                 double z2 = coord[dim_space*N2+2];
1536                 double z3 = coord[dim_space*N3+2];
1537                 double z4 = coord[dim_space*N4+2];
1538
1539                 xbarycenter1 = (x1 + x2 + x3 + x4)/4.0;
1540                 xbarycenter2 = (y1 + y2 + y3 + y4)/4.0;
1541                 xbarycenter3 = (z1 + z2 + z3 + z4)/4.0;
1542                 barycenter[3*index] = xbarycenter1 ;
1543                 barycenter[3*index+1] = xbarycenter2 ;
1544                 barycenter[3*index+2] = xbarycenter3 ;
1545                 index++;
1546               }
1547             break;
1548           }
1549         case MED_PYRA5 : case MED_PYRA13 :
1550           {
1551             for (int pyra=0;pyra<nb_entity_type;pyra++)
1552               {
1553                 int pyra_index = (type%100)*pyra;
1554
1555                 int N1 = global_connectivity[pyra_index]-1;
1556                 int N2 = global_connectivity[pyra_index+1]-1;
1557                 int N3 = global_connectivity[pyra_index+2]-1;
1558                 int N4 = global_connectivity[pyra_index+3]-1;
1559                 int N5 = global_connectivity[pyra_index+4]-1;
1560
1561                 double x1 = coord[dim_space*N1];
1562                 double x2 = coord[dim_space*N2];
1563                 double x3 = coord[dim_space*N3];
1564                 double x4 = coord[dim_space*N4];
1565                 double x5 = coord[dim_space*N5];
1566
1567                 double y1 = coord[dim_space*N1+1];
1568                 double y2 = coord[dim_space*N2+1];
1569                 double y3 = coord[dim_space*N3+1];
1570                 double y4 = coord[dim_space*N4+1];
1571                 double y5 = coord[dim_space*N5+1];
1572
1573                 double z1 = coord[dim_space*N1+2];
1574                 double z2 = coord[dim_space*N2+2];
1575                 double z3 = coord[dim_space*N3+2];
1576                 double z4 = coord[dim_space*N4+2];
1577                 double z5 = coord[dim_space*N5+2];
1578
1579                 xbarycenter1 = (x1 + x2 + x3 + x4 + x5)/5.0;
1580                 xbarycenter2 = (y1 + y2 + y3 + y4 + y5)/5.0;
1581                 xbarycenter3 = (z1 + z2 + z3 + z4 + z5)/5.0;
1582                 barycenter[3*index] = xbarycenter1 ;
1583                 barycenter[3*index+1] = xbarycenter2 ;
1584                 barycenter[3*index+2] = xbarycenter3 ;
1585                 index++;
1586               }
1587             break;
1588           }
1589         case MED_PENTA6 : case MED_PENTA15 :
1590           {
1591             for (int penta=0;penta<nb_entity_type;penta++)
1592               {
1593                 int penta_index = (type%100)*penta;
1594
1595                 int N1 = global_connectivity[penta_index]-1;
1596                 int N2 = global_connectivity[penta_index+1]-1;
1597                 int N3 = global_connectivity[penta_index+2]-1;
1598                 int N4 = global_connectivity[penta_index+3]-1;
1599                 int N5 = global_connectivity[penta_index+4]-1;
1600                 int N6 = global_connectivity[penta_index+5]-1;
1601
1602                 double x1 = coord[dim_space*N1];
1603                 double x2 = coord[dim_space*N2];
1604                 double x3 = coord[dim_space*N3];
1605                 double x4 = coord[dim_space*N4];
1606                 double x5 = coord[dim_space*N5];
1607                 double x6 = coord[dim_space*N6];
1608
1609                 double y1 = coord[dim_space*N1+1];
1610                 double y2 = coord[dim_space*N2+1];
1611                 double y3 = coord[dim_space*N3+1];
1612                 double y4 = coord[dim_space*N4+1];
1613                 double y5 = coord[dim_space*N5+1];
1614                 double y6 = coord[dim_space*N6+1];
1615
1616                 double z1 = coord[dim_space*N1+2];
1617                 double z2 = coord[dim_space*N2+2];
1618                 double z3 = coord[dim_space*N3+2];
1619                 double z4 = coord[dim_space*N4+2];
1620                 double z5 = coord[dim_space*N5+2];
1621                 double z6 = coord[dim_space*N6+2];
1622
1623                 xbarycenter1 = (x1 + x2 + x3 + x4 + x5 + x6)/6.0;
1624                 xbarycenter2 = (y1 + y2 + y3 + y4 + y5 + y6)/6.0;
1625                 xbarycenter3 = (z1 + z2 + z3 + z4 + z5 + z6)/6.0;
1626                 barycenter[3*index] = xbarycenter1 ;
1627                 barycenter[3*index+1] = xbarycenter2 ;
1628                 barycenter[3*index+2] = xbarycenter3 ;
1629                 index++;
1630               }
1631             break;
1632           }
1633         case MED_HEXA8 : case MED_HEXA20 :
1634           {
1635             for (int hexa=0;hexa<nb_entity_type;hexa++)
1636               {
1637                 int hexa_index = (type%100)*hexa;
1638
1639                 int N1 = global_connectivity[hexa_index]-1;
1640                 int N2 = global_connectivity[hexa_index+1]-1;
1641                 int N3 = global_connectivity[hexa_index+2]-1;
1642                 int N4 = global_connectivity[hexa_index+3]-1;
1643                 int N5 = global_connectivity[hexa_index+4]-1;
1644                 int N6 = global_connectivity[hexa_index+5]-1;
1645                 int N7 = global_connectivity[hexa_index+6]-1;
1646                 int N8 = global_connectivity[hexa_index+7]-1;
1647
1648                 double x1 = coord[dim_space*N1];
1649                 double x2 = coord[dim_space*N2];
1650                 double x3 = coord[dim_space*N3];
1651                 double x4 = coord[dim_space*N4];
1652                 double x5 = coord[dim_space*N5];
1653                 double x6 = coord[dim_space*N6];
1654                 double x7 = coord[dim_space*N7];
1655                 double x8 = coord[dim_space*N8];
1656
1657                 double y1 = coord[dim_space*N1+1];
1658                 double y2 = coord[dim_space*N2+1];
1659                 double y3 = coord[dim_space*N3+1];
1660                 double y4 = coord[dim_space*N4+1];
1661                 double y5 = coord[dim_space*N5+1];
1662                 double y6 = coord[dim_space*N6+1];
1663                 double y7 = coord[dim_space*N7+1];
1664                 double y8 = coord[dim_space*N8+1];
1665
1666                 double z1 = coord[dim_space*N1+2];
1667                 double z2 = coord[dim_space*N2+2];
1668                 double z3 = coord[dim_space*N3+2];
1669                 double z4 = coord[dim_space*N4+2];
1670                 double z5 = coord[dim_space*N5+2];
1671                 double z6 = coord[dim_space*N6+2];
1672                 double z7 = coord[dim_space*N7+2];
1673                 double z8 = coord[dim_space*N8+2];
1674
1675                 xbarycenter1 = (x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8)/8.0;
1676                 xbarycenter2 = (y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8)/8.0;
1677                 xbarycenter3 = (z1 + z2 + z3 + z4 + z5 + z6 + z7 + z8)/8.0;
1678
1679                 barycenter[3*index] = xbarycenter1 ;
1680                 barycenter[3*index+1] = xbarycenter2 ;
1681                 barycenter[3*index+2] = xbarycenter3 ;
1682
1683                 index++;
1684               }
1685             break;
1686           }
1687         case MED_TRIA3 : case MED_TRIA6 :
1688           {
1689             for (int tria=0;tria<nb_entity_type;tria++)
1690               {
1691                 int tria_index = (type%100)*tria;
1692
1693                 int N1 = global_connectivity[tria_index]-1;
1694                 int N2 = global_connectivity[tria_index+1]-1;
1695                 int N3 = global_connectivity[tria_index+2]-1;
1696
1697                 double x1 = coord[dim_space*N1];
1698                 double x2 = coord[dim_space*N2];
1699                 double x3 = coord[dim_space*N3];
1700
1701                 double y1 = coord[dim_space*N1+1];
1702                 double y2 = coord[dim_space*N2+1];
1703                 double y3 = coord[dim_space*N3+1];
1704
1705                 xbarycenter1 = (x1 + x2 + x3)/3.0;
1706                 xbarycenter2 = (y1 + y2 + y3)/3.0;
1707
1708                 if (dim_space==2)
1709                   {
1710                     barycenter[2*index] = xbarycenter1 ;
1711                     barycenter[2*index+1] = xbarycenter2 ;
1712                   }
1713                 else
1714                   {
1715                     double z1 =
1716                       coord[dim_space*N1+2];
1717                     double z2 =
1718                       coord[dim_space*N2+2];
1719                     double z3 =
1720                       coord[dim_space*N3+2];
1721
1722                     xbarycenter3 = (z1 + z2 + z3)/3.0;
1723
1724                     barycenter[3*index] = xbarycenter1 ;
1725                     barycenter[3*index+1] = xbarycenter2 ;
1726                     barycenter[3*index+2] = xbarycenter3 ;
1727                   }
1728
1729                 index++;
1730               }
1731             break;
1732           }
1733         case MED_QUAD4 : case MED_QUAD8 :
1734           {
1735             for (int quad=0;quad<nb_entity_type;quad++)
1736               {
1737                 int quad_index = (type%100)*quad;
1738
1739                 int N1 = global_connectivity[quad_index]-1;
1740                 int N2 = global_connectivity[quad_index+1]-1;
1741                 int N3 = global_connectivity[quad_index+2]-1;
1742                 int N4 = global_connectivity[quad_index+3]-1;
1743
1744                 double x1 = coord[dim_space*N1];
1745                 double x2 = coord[dim_space*N2];
1746                 double x3 = coord[dim_space*N3];
1747                 double x4 = coord[dim_space*N4];
1748
1749                 double y1 = coord[dim_space*N1+1];
1750                 double y2 = coord[dim_space*N2+1];
1751                 double y3 = coord[dim_space*N3+1];
1752                 double y4 = coord[dim_space*N4+1];
1753
1754                 xbarycenter1 = (x1 + x2 + x3 + x4)/4.0;
1755                 xbarycenter2 = (y1 + y2 + y3 + y4)/4.0;
1756
1757                 if (dim_space==2)
1758                   {
1759                     barycenter[2*index] = xbarycenter1 ;
1760                     barycenter[2*index+1] = xbarycenter2 ;
1761                   }
1762                 else
1763                   {
1764                     double z1 = coord[dim_space*N1+2];
1765                     double z2 = coord[dim_space*N2+2];
1766                     double z3 = coord[dim_space*N3+2];
1767                     double z4 = coord[dim_space*N4+2];
1768
1769                     xbarycenter3 = (z1 + z2 + z3 + z4)/4.0;
1770
1771                     barycenter[3*index] = xbarycenter1 ;
1772                     barycenter[3*index+1] = xbarycenter2 ;
1773                     barycenter[3*index+2] = xbarycenter3 ;
1774                   }
1775                 index++;
1776               }
1777             break;
1778           }
1779         case MED_SEG2 : case MED_SEG3 :
1780           {
1781             for (int edge=0;edge<nb_entity_type;edge++)
1782               {
1783                 int edge_index = (type%100)*edge;
1784
1785                 int N1 = global_connectivity[edge_index]-1;
1786                 int N2 = global_connectivity[edge_index+1]-1;
1787
1788                 double x1 = coord[dim_space*N1];
1789                 double x2 = coord[dim_space*N2];
1790
1791                 double y1 = coord[dim_space*N1+1];
1792                 double y2 = coord[dim_space*N2+1];
1793
1794                 xbarycenter1 = (x1 + x2)/2.0;
1795                 xbarycenter2 = (y1 + y2)/2.0;
1796
1797                 if (dim_space==2)
1798                   {
1799                     barycenter[2*index] = xbarycenter1 ;
1800                     barycenter[2*index+1] = xbarycenter2 ;
1801                   }
1802                 else
1803                   {
1804                     double z1 = coord[dim_space*N1+2];
1805                     double z2 = coord[dim_space*N2+2];
1806
1807                     xbarycenter3 = (z1 + z2)/2.0;
1808
1809                     barycenter[3*index] = xbarycenter1 ;
1810                     barycenter[3*index+1] = xbarycenter2 ;
1811                     barycenter[3*index+2] = xbarycenter3 ;
1812                   }
1813                 index++;
1814               }
1815             break;
1816           }
1817         default :
1818           throw MEDEXCEPTION("MESH::getBarycenter(SUPPORT*) Bad Support to get a barycenter on it (in fact unknown type) !");
1819           break;
1820         }
1821     }
1822
1823   Barycenter->setValue(MED_FULL_INTERLACE,barycenter);
1824
1825   return Barycenter;
1826 }