src/MEDMEM/MEDMEM_Mesh.cxx

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