src/MEDMEM/MEDMEM_Mesh.cxx

   1 //  MED MEDMEM : MED files in memory
   2 //
   3 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  21 //
  22 //
  23 //
  24 //  File   : MEDMEM_Mesh.cxx
  25 //  Module : MED
  26
  27 using namespace std;
  28 /*
  29  File Mesh.cxx
  30  $Header$
  31 */
  32
  33 #include <math.h>
  34
  35 #include <list>
  36 #include <map>
  37
  38 #include "MEDMEM_DriversDef.hxx"
  39 #include "MEDMEM_Field.hxx"
  40 #include "MEDMEM_Mesh.hxx"
  41
  42 #include "MEDMEM_Support.hxx"
  43 #include "MEDMEM_Family.hxx"
  44 #include "MEDMEM_Group.hxx"
  45 #include "MEDMEM_Coordinate.hxx"
  46 #include "MEDMEM_Connectivity.hxx"
  47 #include "MEDMEM_CellModel.hxx"
  48 #include "MEDMEM_Grid.hxx"
  49
  50 //update Families with content list
  51 //int family_count(int family_number, int count, int * entities_number, int * entities_list) ;
  52
  53 // ------- Drivers Management Part
  54
  55 // MESH::INSTANCE_DE<MED_MESH_RDONLY_DRIVER> MESH::inst_med_rdonly ;
  56 //const MESH::INSTANCE * const MESH::instances[] = { &MESH::inst_med_rdonly , &MESH::inst_med_rdwr } ;
  57
  58 // Add a similar line for your personnal driver (step 3)
  59 MESH::INSTANCE_DE<MED_MESH_RDWR_DRIVER>   MESH::inst_med   ;
  60 MESH::INSTANCE_DE<GIBI_MESH_RDWR_DRIVER>  MESH::inst_gibi ;
  61 // Add your own driver in the driver list       (step 4)
  62 // Note the list must be coherent with the driver type list defined in MEDMEM_DRIVER.hxx.
  63 const MESH::INSTANCE * const MESH::instances[] =   {  &MESH::inst_med, &MESH::inst_gibi } ;
  64
  65 /*! Add a MESH driver of type (MED_DRIVER, ....) associated with file <fileName>. The meshname used in the file
  66     is  <driverName>. addDriver returns an int handler. */
  67 int MESH::addDriver(driverTypes driverType,
  68                     const string & fileName/*="Default File Name.med"*/,const string & driverName/*="Default Mesh Name"*/) {
  69
  70   const char * LOC = "MESH::addDriver(driverTypes driverType, const string & fileName=\"Default File Name.med\",const string & driverName=\"Default Mesh Name\") : ";
  71
  72   GENDRIVER * driver;
  73   int current;
  74
  75   BEGIN_OF(LOC);
  76
  77   driver = instances[driverType]->run(fileName, this) ;
  78   _drivers.push_back(driver);
  79   current = _drivers.size()-1;
  80
  81   _drivers[current]->setMeshName(driverName);
  82   return current;
  83
  84   END_OF(LOC);
  85
  86 }
  87
  88 /*! Add an existing MESH driver. */
  89 int  MESH::addDriver(MED_MESH_DRIVER & driver) {
  90   const char * LOC = "MESH::addDriver(GENDRIVER &) : ";
  91   BEGIN_OF(LOC);
  92
  93   // A faire : Vérifier que le driver est de type MESH.
  94   GENDRIVER * newDriver = driver.copy() ;
  95
  96   _drivers.push_back(newDriver);
  97   return _drivers.size()-1;
  98
  99   END_OF(LOC);
 100 }
 101
 102 /*! Remove an existing MESH driver. */
 103 void MESH::rmDriver (int index/*=0*/) {
 104   const char * LOC = "MESH::rmDriver (int index=0): ";
 105   BEGIN_OF(LOC);
 106
 107   if ( _drivers[index] ) {
 108     //_drivers.erase(&_drivers[index]);
 109     // why not ????
 110     MESSAGE ("detruire");
 111   }
 112   else
 113     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
 114                                      << "The index given is invalid, index must be between  0 and  |"
 115                                      << _drivers.size()
 116                                      )
 117                           );
 118
 119   END_OF(LOC);
 120
 121 };
 122
 123 // ------ End of Drivers Management Part
 124
 125
 126 void MESH::init() {
 127
 128   const char * LOC = "MESH::init(): ";
 129
 130   BEGIN_OF(LOC);
 131
 132   string        _name = "NOT DEFINED"; // A POSITIONNER EN FCT DES IOS ?
 133
 134   COORDINATE   * _coordinate   = (COORDINATE   *) NULL;
 135   CONNECTIVITY * _connectivity = (CONNECTIVITY *) NULL;
 136
 137   _spaceDimension        =          MED_INVALID; // 0 ?!?
 138   _meshDimension         =          MED_INVALID;
 139   _numberOfNodes         =          MED_INVALID;
 140
 141   _numberOfNodesFamilies =          0; // MED_INVALID ?!?
 142   _numberOfCellsFamilies =          0;
 143   _numberOfFacesFamilies =          0;
 144   _numberOfEdgesFamilies =          0;
 145
 146   _numberOfNodesGroups   =          0;  // MED_INVALID ?!?
 147   _numberOfCellsGroups   =          0;
 148   _numberOfFacesGroups   =          0;
 149   _numberOfEdgesGroups   =          0;
 150
 151   _isAGrid = false;
 152 };
 153
 154 /*! Create an empty MESH. */
 155 MESH::MESH():_coordinate(NULL),_connectivity(NULL) {
 156   init();
 157 };
 158
 159 MESH::MESH(MESH &m)
 160 {
 161   _name=m._name;
 162   _isAGrid = m._isAGrid;
 163
 164   if (m._coordinate != NULL)
 165     _coordinate = new COORDINATE(* m._coordinate);
 166   else
 167     _coordinate = (COORDINATE *) NULL;
 168
 169   if (m._connectivity != NULL)
 170     _connectivity = new CONNECTIVITY(* m._connectivity);
 171   else
 172     _connectivity = (CONNECTIVITY *) NULL;
 173
 174   _spaceDimension = m._spaceDimension;
 175   _meshDimension = m._meshDimension;
 176   _numberOfNodes = m._numberOfNodes;
 177
 178   _numberOfNodesFamilies = m._numberOfNodesFamilies;
 179   _familyNode = m._familyNode;
 180   for (int i=0; i<m._numberOfNodesFamilies; i++)
 181     {
 182       _familyNode[i] = new FAMILY(* m._familyNode[i]);
 183       _familyNode[i]->setMesh(this);
 184     }
 185
 186   _numberOfCellsFamilies = m._numberOfCellsFamilies;
 187   _familyCell = m._familyCell;
 188   for (int i=0; i<m._numberOfCellsFamilies; i++)
 189     {
 190       _familyCell[i] = new FAMILY(* m._familyCell[i]);
 191       _familyCell[i]->setMesh(this);
 192     }
 193
 194   _numberOfFacesFamilies = m._numberOfFacesFamilies;
 195   _familyFace = m._familyFace;
 196   for (int i=0; i<m._numberOfFacesFamilies; i++)
 197     {
 198       _familyFace[i] = new FAMILY(* m._familyFace[i]);
 199       _familyFace[i]->setMesh(this);
 200     }
 201
 202   _numberOfEdgesFamilies = m._numberOfEdgesFamilies;
 203   _familyEdge = m._familyEdge;
 204   for (int i=0; i<m._numberOfEdgesFamilies; i++)
 205     {
 206       _familyEdge[i] = new FAMILY(* m._familyEdge[i]);
 207       _familyEdge[i]->setMesh(this);
 208     }
 209
 210   _numberOfNodesGroups = m._numberOfNodesGroups;
 211   _groupNode = m._groupNode;
 212   for (int i=0; i<m._numberOfNodesGroups; i++)
 213     {
 214       _groupNode[i] = new GROUP(* m._groupNode[i]);
 215       _groupNode[i]->setMesh(this);
 216     }
 217
 218   _numberOfCellsGroups = m._numberOfCellsGroups;
 219   _groupCell = m._groupCell;
 220   for (int i=0; i<m._numberOfCellsGroups; i++)
 221     {
 222       _groupCell[i] = new GROUP(* m._groupCell[i]);
 223       _groupCell[i]->setMesh(this);
 224     }
 225
 226   _numberOfFacesGroups = m._numberOfFacesGroups;
 227   _groupFace = m._groupFace;
 228   for (int i=0; i<m._numberOfFacesGroups; i++)
 229     {
 230       _groupFace[i] = new GROUP(* m._groupFace[i]);
 231       _groupFace[i]->setMesh(this);
 232     }
 233
 234   _numberOfEdgesGroups = m._numberOfEdgesGroups;
 235   _groupEdge = m._groupEdge;
 236   for (int i=0; i<m._numberOfEdgesGroups; i++)
 237     {
 238       _groupEdge[i] = new GROUP(* m._groupEdge[i]);
 239       _groupEdge[i]->setMesh(this);
 240     }
 241
 242   //_drivers = m._drivers;  //Recopie des drivers?
 243 }
 244
 245 MESH::~MESH() {
 246
 247   MESSAGE("MESH::~MESH() : Destroying the Mesh");
 248   if (_coordinate != ((COORDINATE *) NULL)) delete _coordinate ;
 249   if (_connectivity != ((CONNECTIVITY *) NULL)) delete _connectivity ;
 250   int size ;
 251   size = _familyNode.size() ;
 252   for (int i=0;i<size;i++)
 253     ;// mpv such destructors call is problemmatic for ALLIANCES algorithms
 254     //delete _familyNode[i] ;
 255   size = _familyCell.size() ;
 256   for (int i=0;i<size;i++)
 257     ;// mpv such destructors call is problemmatic for ALLIANCES algorithms
 258     //delete _familyCell[i] ;
 259   size = _familyFace.size() ;
 260   for (int i=0;i<size;i++)
 261     ;// mpv such destructors call is problemmatic for ALLIANCES algorithms
 262     //delete _familyFace[i] ;
 263   size = _familyEdge.size() ;
 264   for (int i=0;i<size;i++)
 265     ;// mpv such destructors call is problemmatic for ALLIANCES algorithms
 266     //delete _familyEdge[i] ;
 267   size = _groupNode.size() ;
 268   for (int i=0;i<size;i++)
 269     ;// mpv such destructors call is problemmatic for ALLIANCES algorithms
 270     //delete _groupNode[i] ;
 271   size = _groupCell.size() ;
 272   for (int i=0;i<size;i++)
 273     ;// mpv such destructors call is problemmatic for ALLIANCES algorithms
 274     //delete _groupCell[i] ;
 275   size = _groupFace.size() ;
 276   for (int i=0;i<size;i++)
 277     ;// mpv such destructors call is problemmatic for ALLIANCES algorithms
 278     //delete _groupFace[i] ;
 279   size = _groupEdge.size() ;
 280   for (int i=0;i<size;i++)
 281     ;// mpv such destructors call is problemmatic for ALLIANCES algorithms
 282     //delete _groupEdge[i] ;
 283
 284   MESSAGE("In this object MESH there is(are) " << _drivers.size() << " driver(s)");
 285
 286   for (int index=0; index < _drivers.size(); index++ )
 287     {
 288       SCRUTE(_drivers[index]);
 289       if ( _drivers[index] != NULL) delete _drivers[index];
 290     }
 291
 292 }
 293
 294 MESH & MESH::operator=(const MESH &m)
 295 {
 296
 297   //  A FAIRE.........
 298
 299   // ATTENTION CET OPERATEUR DE RECOPIE EST DANGEREUX POUR LES
 300   // POINTEURS : ex : nodal_connectivity ???? EXPRES ????
 301 //    _drivers = m._drivers;
 302
 303 //    space_dimension=m.space_dimension;
 304 //    mesh_dimension=m.mesh_dimension;
 305
 306 //    nodes_count=m.nodes_count;
 307 //    coordinates=m.coordinates;
 308 //    coordinates_name=m.coordinates_name ;
 309 //    coordinates_unit=m.coordinates_unit ;
 310 //    nodes_numbers=m.nodes_numbers ;
 311
 312 //    cells_types_count=m.cells_types_count;
 313 //    cells_type=m.cells_type;
 314 //    cells_count=m.cells_count;
 315 //    nodal_connectivity=m.nodal_connectivity;
 316
 317 //    nodes_families_count=m.nodes_families_count;
 318 //    nodes_Families=m.nodes_Families;
 319
 320 //    cells_families_count=m.cells_families_count;
 321 //    cells_Families=m.cells_Families;
 322
 323 //    maximum_cell_number_by_node = m.maximum_cell_number_by_node;
 324 //    if (maximum_cell_number_by_node > 0)
 325 //      {
 326 //        reverse_nodal_connectivity = m.reverse_nodal_connectivity;
 327 //        reverse_nodal_connectivity_index = m.reverse_nodal_connectivity_index ;
 328 //      }
 329   return *this;
 330 }
 331
 332 /*! Create a MESH object using a MESH driver of type (MED_DRIVER, ....) associated with file <fileName>.
 333   The meshname <driverName> must already exists in the file.*/
 334 MESH::MESH(driverTypes driverType, const string &  fileName/*=""*/, const string &  driverName/*=""*/) {
 335   const char * LOC ="MESH::MESH(driverTypes driverType, const string &  fileName="", const string &  driverName="") : ";
 336
 337   int current;
 338
 339   BEGIN_OF(LOC);
 340
 341   init();
 342
 343   MED_MESH_RDONLY_DRIVER myDriver(fileName,this) ;
 344   myDriver.setMeshName(driverName);
 345   current = addDriver(myDriver);
 346 //   current = addDriver(driverType,fileName,driverName);
 347 //   switch(_drivers[current]->getAccessMode() ) {
 348 //   case MED_WRONLY : {
 349 //     MESSAGE("MESH::MESH(driverTypes driverType, .....) : driverType must not be a MED_WRONLY accessMode");
 350 //     rmDriver(current);
 351 //     break;}
 352 //   default : {
 353 //   }
 354 //   }
 355   _drivers[current]->open();
 356   _drivers[current]->read();
 357   _drivers[current]->close();
 358
 359   if (_isAGrid)
 360     ((GRID *) this)->fillMeshAfterRead();
 361
 362   END_OF(LOC);
 363 };
 364
 365
 366 //  Node MESH::Donne_Barycentre(const Maille &m) const
 367 //  {
 368 //    Node N=node[m[0]];
 369 //    for (int i=1;i<m.donne_nbr_sommets();i++) N=N+node[m[i]];
 370 //    N=N*(1.0/m.donne_nbr_sommets());
 371 //    return N;
 372 //  }
 373
 374 ostream & operator<<(ostream &os, MESH &myMesh)
 375 {
 376   int spacedimension = myMesh.getSpaceDimension();
 377   int meshdimension  = myMesh.getMeshDimension();
 378   int numberofnodes  = myMesh.getNumberOfNodes();
 379
 380   os << "Space Dimension : " << spacedimension << endl << endl;
 381
 382   os << "Mesh Dimension : " << meshdimension << endl << endl;
 383
 384   const double * coordinates = myMesh.getCoordinates(MED_FULL_INTERLACE);
 385   os << "SHOW NODES COORDINATES : " << endl;
 386
 387   os << "Name :" << endl;
 388   const string * coordinatesnames = myMesh.getCoordinatesNames();
 389   for(int i=0; i<spacedimension ; i++)
 390     {
 391       os << " - " << coordinatesnames[i] << endl;
 392     }
 393   os << "Unit :" << endl;
 394   const string * coordinatesunits = myMesh.getCoordinatesUnits();
 395   for(int i=0; i<spacedimension ; i++)
 396     {
 397       os << " - " << coordinatesunits[i] << endl;
 398     }
 399   for(int i=0; i<numberofnodes ; i++)
 400     {
 401       os << "Nodes " << i+1 << " : ";
 402       for (int j=0; j<spacedimension ; j++)
 403         os << coordinates[i*spacedimension+j] << " ";
 404       os << endl;
 405     }
 406
 407   const CONNECTIVITY * myConnectivity = myMesh.getConnectivityptr();
 408   if (!myConnectivity->existConnectivity(MED_NODAL,MED_CELL) && myConnectivity->existConnectivity(MED_DESCENDING,MED_CELL))
 409     {
 410       os << endl << "SHOW CONNECTIVITY (DESCENDING) :" << endl;
 411       int numberofelements;
 412       const int * connectivity;
 413       const int * connectivity_index;
 414       myMesh.calculateConnectivity(MED_FULL_INTERLACE,MED_DESCENDING,MED_CELL);
 415       try {
 416         numberofelements = myMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
 417         connectivity =  myMesh.getConnectivity(MED_FULL_INTERLACE,MED_DESCENDING,MED_CELL,MED_ALL_ELEMENTS);
 418         connectivity_index =  myMesh.getConnectivityIndex(MED_DESCENDING,MED_CELL);
 419       }
 420       catch (MEDEXCEPTION m) {
 421         os << m.what() << endl;
 422         exit (-1);
 423       }
 424       for (int j=0;j<numberofelements;j++) {
 425         os << "Element "<<j+1<<" : ";
 426         for (int k=connectivity_index[j];k<connectivity_index[j+1];k++)
 427           os << connectivity[k-1]<<" ";
 428         os << endl;
 429       }
 430     }
 431   else
 432     {
 433       int numberoftypes = myMesh.getNumberOfTypes(MED_CELL);
 434       const medGeometryElement  * types = myMesh.getTypes(MED_CELL);
 435       os << endl << "SHOW CONNECTIVITY (NODAL) :" << endl;
 436       for (int i=0; i<numberoftypes; i++) {
 437         os << "For type " << types[i] << " : " << endl;
 438         int numberofelements = myMesh.getNumberOfElements(MED_CELL,types[i]);
 439         const int * connectivity =  myMesh.getConnectivity(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,types[i]);
 440         int numberofnodespercell = types[i]%100;
 441         for (int j=0;j<numberofelements;j++){
 442           os << "Element "<< j+1 <<" : ";
 443           for (int k=0;k<numberofnodespercell;k++)
 444             os << connectivity[j*numberofnodespercell+k]<<" ";
 445           os << endl;
 446         }
 447       }
 448     }
 449
 450
 451   medEntityMesh entity;
 452   os << endl << "SHOW FAMILIES :" << endl << endl;
 453   for (int k=1; k<=4; k++)
 454     {
 455       if (k==1) entity = MED_NODE;
 456       if (k==2) entity = MED_CELL;
 457       if (k==3) entity = MED_FACE;
 458       if (k==4) entity = MED_EDGE;
 459       int numberoffamilies = myMesh.getNumberOfFamilies(entity);
 460       using namespace MED_FR;
 461       os << "NumberOfFamilies on "<< entNames[(MED_FR::med_entite_maillage) entity] <<" : "<<numberoffamilies<<endl;
 462       using namespace MED_EN;
 463       for (int i=1; i<numberoffamilies+1;i++)
 464         {
 465           os << * myMesh.getFamily(entity,i) << endl;
 466         }
 467     }
 468
 469   os << endl << "SHOW GROUPS :" << endl << endl;
 470   for (int k=1; k<=4; k++)
 471     {
 472       if (k==1) entity = MED_NODE;
 473       if (k==2) entity = MED_CELL;
 474       if (k==3) entity = MED_FACE;
 475       if (k==4) entity = MED_EDGE;
 476       int numberofgroups = myMesh.getNumberOfGroups(entity);
 477       using namespace MED_FR;
 478       os << "NumberOfGroups on "<< entNames[(MED_FR::med_entite_maillage) entity] <<" : "<<numberofgroups<<endl;
 479       using namespace MED_EN;
 480       for (int i=1; i<numberofgroups+1;i++)
 481         {
 482           os << * myMesh.getGroup(entity,i) << endl;
 483         }
 484     }
 485
 486   return os;
 487 }
 488
 489 /*!
 490   Get global number of element which have same connectivity than connectivity argument.
 491
 492   It do not take care of connectivity order (3,4,7,10 is same as 7,3,10,4).
 493
 494   Return -1 if not found.
 495 */
 496 int MESH::getElementNumber(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type, int * connectivity) const
 497 {
 498   const char* LOC="MESH::getElementNumber " ;
 499   BEGIN_OF(LOC) ;
 500
 501   int numberOfValue ; // size of connectivity array
 502   CELLMODEL myType(Type) ;
 503   if (ConnectivityType==MED_DESCENDING)
 504     numberOfValue = myType.getNumberOfConstituents(1) ; // faces (3D) or edges (2D)
 505   else
 506     numberOfValue = myType.getNumberOfNodes() ; // nodes
 507
 508   const int * myReverseConnectivityValue = getReverseConnectivity(ConnectivityType,Entity) ;
 509   const int * myReverseConnectivityIndex = getReverseConnectivityIndex(ConnectivityType,Entity) ;
 510
 511   // First node or face/edge
 512   int indexBegin = myReverseConnectivityIndex[connectivity[0]-1] ;
 513   int indexEnd = myReverseConnectivityIndex[connectivity[0]] ;
 514
 515   // list to put cells numbers
 516   list<int> cellsList ;
 517   list<int>::iterator itList ;
 518   for (int i=indexBegin; i<indexEnd; i++)
 519     cellsList.push_back(myReverseConnectivityValue[i-1]) ;
 520
 521   // Foreach node or face/edge in cell, we search cells which are in common.
 522   // At the end, we must have only one !
 523
 524   for (int i=1; i<numberOfValue; i++) {
 525     int connectivity_i = connectivity[i] ;
 526     for (itList=cellsList.begin();itList!=cellsList.end();itList++) {
 527       bool find = false ;
 528       for (int j=myReverseConnectivityIndex[connectivity_i-1]; j<myReverseConnectivityIndex[connectivity_i]; j++) {
 529         if ((*itList)==myReverseConnectivityValue[j-1]) {
 530           find=true ;
 531           break ;
 532         }
 533       }
 534       if (!find) {
 535         itList=cellsList.erase(itList);
 536         itList--; // well : rigth if itList = cellsList.begin() ??
 537       }
 538     }
 539   }
 540
 541   if (cellsList.size()>1) // we have more than one cell
 542     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Found more than one element !")) ;
 543
 544   if (cellsList.size()==0)
 545     return -1;
 546
 547   END_OF(LOC);
 548
 549   return cellsList.front() ;
 550 }
 551
 552
 553 /*!
 554   Return a support which reference all elements on the boundary of mesh.
 555
 556   For instance, we could get only face in 3D and edge in 2D.
 557 */
 558 SUPPORT * MESH::getBoundaryElements(medEntityMesh Entity)
 559   throw (MEDEXCEPTION)
 560 {
 561   const char * LOC = "MESH::getBoundaryElements : " ;
 562   BEGIN_OF(LOC) ;
 563   // some test :
 564   // actually we could only get face (in 3D) and edge (in 2D)
 565   if (_spaceDimension == 3)
 566     if (Entity != MED_FACE)
 567       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<Entity<<" !"));
 568   if (_spaceDimension == 2)
 569     if (Entity != MED_EDGE)
 570       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<Entity<<" !"));
 571
 572   // well, all rigth !
 573   SUPPORT * mySupport = new SUPPORT(this,"Boundary",Entity);
 574   //mySupport.setDescription("boundary of type ...");
 575   mySupport->setAll(false);
 576
 577
 578   const int * myConnectivityValue = getReverseConnectivity(MED_DESCENDING) ;
 579   const int * myConnectivityIndex = getReverseConnectivityIndex(MED_DESCENDING) ;
 580   int numberOf = getNumberOfElements(Entity,MED_ALL_ELEMENTS) ;
 581   list<int> myElementsList ;
 582   int size = 0 ;
 583   for (int i=0 ; i<numberOf; i++)
 584     if (myConnectivityValue[myConnectivityIndex[i]] == 0) {
 585       myElementsList.push_back(i+1) ;
 586       size++ ;
 587     }
 588   // Well, we must know how many geometric type we have found
 589   int * myListArray = new int[size] ;
 590   int id = 0 ;
 591   list<int>::iterator myElementsListIt ;
 592   for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
 593     myListArray[id]=(*myElementsListIt) ;
 594     id ++ ;
 595   }
 596
 597   int numberOfGeometricType ;
 598   medGeometryElement* geometricType ;
 599   int * numberOfGaussPoint ;
 600   int * geometricTypeNumber ;
 601   int * numberOfElements ;
 602   //MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY() ;
 603   int * mySkyLineArrayIndex ;
 604
 605   int numberOfType = getNumberOfTypes(Entity) ;
 606   if (numberOfType == 1) { // wonderfull : it's easy !
 607     numberOfGeometricType = 1 ;
 608     geometricType = new medGeometryElement[1] ;
 609     const medGeometryElement *  allType = getTypes(Entity);
 610     geometricType[0] = allType[0] ;
 611     numberOfGaussPoint = new int[1] ;
 612     numberOfGaussPoint[0] = 1 ;
 613     geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
 614     geometricTypeNumber[0] = 0 ;
 615     numberOfElements = new int[1] ;
 616     numberOfElements[0] = size ;
 617     mySkyLineArrayIndex = new int[2] ;
 618     mySkyLineArrayIndex[0]=1 ;
 619     mySkyLineArrayIndex[1]=1+size ;
 620   }
 621   else {// hemmm
 622     map<medGeometryElement,int> theType ;
 623     for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
 624       medGeometryElement myType = getElementType(Entity,*myElementsListIt) ;
 625       if (theType.find(myType) != theType.end() )
 626         theType[myType]+=1 ;
 627       else
 628         theType[myType]=1 ;
 629     }
 630     numberOfGeometricType = theType.size() ;
 631     geometricType = new medGeometryElement[numberOfGeometricType] ;
 632     const medGeometryElement *  allType = getTypes(Entity);
 633     numberOfGaussPoint = new int[numberOfGeometricType] ;
 634     geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
 635     numberOfElements = new int[numberOfGeometricType] ;
 636     mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
 637     int index = 0 ;
 638     mySkyLineArrayIndex[0]=1 ;
 639     map<medGeometryElement,int>::iterator theTypeIt ;
 640     for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
 641       geometricType[index] = (*theTypeIt).first ;
 642       numberOfGaussPoint[index] = 1 ;
 643       geometricTypeNumber[index] = 0 ;
 644       numberOfElements[index] = (*theTypeIt).second ;
 645       mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfElements[index] ;
 646       index++ ;
 647     }
 648   }
 649   //mySkyLineArray->setMEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
 650   MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
 651
 652   mySupport->setNumberOfGeometricType(numberOfGeometricType) ;
 653   mySupport->setGeometricType(geometricType) ;
 654   mySupport->setNumberOfGaussPoint(numberOfGaussPoint) ;
 655   mySupport->setNumberOfElements(numberOfElements) ;
 656   mySupport->setTotalNumberOfElements(size) ;
 657   // mySupport->setGeometricTypeNumber(geometricTypeNumber) ;
 658   mySupport->setNumber(mySkyLineArray) ;
 659
 660   delete[] numberOfElements;
 661   delete[] geometricTypeNumber;
 662   delete[] numberOfGaussPoint;
 663   delete[] geometricType;
 664   delete[] mySkyLineArrayIndex;
 665   delete[] myListArray;
 666 //   delete mySkyLineArray;
 667
 668   END_OF(LOC) ;
 669   return mySupport ;
 670 }
 671
 672 FIELD<double>* MESH::getVolume(const SUPPORT * Support) const throw (MEDEXCEPTION)
 673 {
 674   const char * LOC = "MESH::getVolume(SUPPORT*) : ";
 675   BEGIN_OF(LOC);
 676
 677   // Support must be on 3D elements
 678
 679   // Make sure that the MESH class is the same as the MESH class attribut
 680   // in the class Support
 681   MESH* myMesh = Support->getMesh();
 682   if (this != myMesh)
 683     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
 684
 685   int dim_space = getSpaceDimension();
 686   medEntityMesh support_entity = Support->getEntity();
 687   bool onAll = Support->isOnAllElements();
 688
 689   int nb_type, length_values;
 690   const medGeometryElement* types;
 691   int nb_entity_type;
 692   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
 693   const int* global_connectivity;
 694
 695 //    if (onAll)
 696 //      {
 697 //        nb_type = myMesh->getNumberOfTypes(support_entity);
 698 //        length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
 699 //        types = getTypes(support_entity);
 700 //      }
 701 //    else
 702     {
 703       nb_type = Support->getNumberOfTypes();
 704       length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
 705       types = Support->getTypes();
 706     }
 707
 708   int index;
 709   FIELD<double>* Volume = new FIELD<double>::FIELD(Support,1);
 710   //  double *volume = new double [length_values];
 711   Volume->setName("VOLUME");
 712   Volume->setDescription("cells volume");
 713   Volume->setComponentName(1,"volume");
 714   Volume->setComponentDescription(1,"desc-comp");
 715
 716   /*  string MEDComponentUnit(MED_TAILLE_PNOM,' ');*/
 717
 718   string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1]+"*"+myMesh->getCoordinatesUnits()[2];
 719
 720   Volume->setMEDComponentUnit(1,MEDComponentUnit);
 721
 722   Volume->setValueType(MED_REEL64);
 723
 724   Volume->setIterationNumber(0);
 725   Volume->setOrderNumber(0);
 726   Volume->setTime(0.0);
 727
 728   //const double *volume = Volume->getValue(MED_FULL_INTERLACE);
 729   MEDARRAY<double> *volume = Volume->getvalue();
 730
 731   index = 1;
 732   const double * coord = getCoordinates(MED_FULL_INTERLACE);
 733
 734   for (int i=0;i<nb_type;i++)
 735     {
 736       medGeometryElement type = types[i] ;
 737       double xvolume;
 738
 739       if (onAll)
 740         {
 741           nb_entity_type = getNumberOfElements(support_entity,type);
 742           global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
 743         }
 744       else
 745         {
 746           // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all"));
 747
 748           nb_entity_type = Support->getNumberOfElements(type);
 749
 750           const int * supp_number = Support->getNumber(type);
 751           const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
 752           const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
 753           int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
 754
 755           for (int k_type = 0; k_type<nb_entity_type; k_type++) {
 756             for (int j_ent = 0; j_ent<(type%100); j_ent++) {
 757               global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
 758             }
 759           }
 760           global_connectivity = global_connectivity_tmp ;
 761         }
 762
 763       switch (type)
 764         {
 765         case MED_TETRA4 : case MED_TETRA10 :
 766           {
 767             for (int tetra=0;tetra<nb_entity_type;tetra++)
 768               {
 769                 int tetra_index = (type%100)*tetra;
 770
 771                 int N1 = global_connectivity[tetra_index]-1;
 772                 int N2 = global_connectivity[tetra_index+1]-1;
 773                 int N3 = global_connectivity[tetra_index+2]-1;
 774                 int N4 = global_connectivity[tetra_index+3]-1;
 775
 776                 double x1 = coord[dim_space*N1];
 777                 double x2 = coord[dim_space*N2];
 778                 double x3 = coord[dim_space*N3];
 779                 double x4 = coord[dim_space*N4];
 780
 781                 double y1 = coord[dim_space*N1+1];
 782                 double y2 = coord[dim_space*N2+1];
 783                 double y3 = coord[dim_space*N3+1];
 784                 double y4 = coord[dim_space*N4+1];
 785
 786                 double z1 = coord[dim_space*N1+2];
 787                 double z2 = coord[dim_space*N2+2];
 788                 double z3 = coord[dim_space*N3+2];
 789                 double z4 = coord[dim_space*N4+2];
 790
 791                 xvolume = ((x3-x1)*((y2-y1)*(z4-z1) - (z2-z1)*(y4-y1)) -
 792                            (x2-x1)*((y3-y1)*(z4-z1) - (z3-z1)*(y4-y1)) +
 793                            (x4-x1)*((y3-y1)*(z2-z1) - (z3-z1)*(y2-y1)))/6.0;
 794
 795                 //volume[index] = xvolume ;
 796                 volume->setIJ(index,1,xvolume) ;
 797                 index++;
 798               }
 799             break;
 800           }
 801         case MED_PYRA5 : case MED_PYRA13 :
 802           {
 803             for (int pyra=0;pyra<nb_entity_type;pyra++)
 804               {
 805                 int pyra_index = (type%100)*pyra;
 806
 807                 int N1 = global_connectivity[pyra_index]-1;
 808                 int N2 = global_connectivity[pyra_index+1]-1;
 809                 int N3 = global_connectivity[pyra_index+2]-1;
 810                 int N4 = global_connectivity[pyra_index+3]-1;
 811                 int N5 = global_connectivity[pyra_index+4]-1;
 812
 813                 double x1 = coord[dim_space*N1];
 814                 double x2 = coord[dim_space*N2];
 815                 double x3 = coord[dim_space*N3];
 816                 double x4 = coord[dim_space*N4];
 817                 double x5 = coord[dim_space*N5];
 818
 819                 double y1 = coord[dim_space*N1+1];
 820                 double y2 = coord[dim_space*N2+1];
 821                 double y3 = coord[dim_space*N3+1];
 822                 double y4 = coord[dim_space*N4+1];
 823                 double y5 = coord[dim_space*N5+1];
 824
 825                 double z1 = coord[dim_space*N1+2];
 826                 double z2 = coord[dim_space*N2+2];
 827                 double z3 = coord[dim_space*N3+2];
 828                 double z4 = coord[dim_space*N4+2];
 829                 double z5 = coord[dim_space*N5+2];
 830
 831                 xvolume = (((x3-x1)*((y2-y1)*(z5-z1) - (z2-z1)*(y5-y1)) -
 832                             (x2-x1)*((y3-y1)*(z5-z1) - (z3-z1)*(y5-y1)) +
 833                             (x5-x1)*((y3-y1)*(z2-z1) - (z3-z1)*(y2-y1))) +
 834                            ((x4-x1)*((y3-y1)*(z5-z1) - (z3-z1)*(y5-y1)) -
 835                             (x3-x1)*((y4-y1)*(z5-z1) - (z4-z1)*(y5-y1)) +
 836                             (x5-x1)*((y4-y1)*(z3-z1) - (z4-z1)*(y3-y1)))
 837                            )/6.0;
 838
 839                 //volume[index] = xvolume ;
 840                 volume->setIJ(index,1,xvolume) ;
 841                 index = index++;
 842               }
 843             break;
 844           }
 845         case MED_PENTA6 : case MED_PENTA15 :
 846           {
 847             for (int penta=0;penta<nb_entity_type;penta++)
 848               {
 849                 int penta_index = (type%100)*penta;
 850
 851                 int N1 = global_connectivity[penta_index]-1;
 852                 int N2 = global_connectivity[penta_index+1]-1;
 853                 int N3 = global_connectivity[penta_index+2]-1;
 854                 int N4 = global_connectivity[penta_index+3]-1;
 855                 int N5 = global_connectivity[penta_index+4]-1;
 856                 int N6 = global_connectivity[penta_index+5]-1;
 857
 858                 double x1 = coord[dim_space*N1];
 859                 double x2 = coord[dim_space*N2];
 860                 double x3 = coord[dim_space*N3];
 861                 double x4 = coord[dim_space*N4];
 862                 double x5 = coord[dim_space*N5];
 863                 double x6 = coord[dim_space*N6];
 864
 865                 double y1 = coord[dim_space*N1+1];
 866                 double y2 = coord[dim_space*N2+1];
 867                 double y3 = coord[dim_space*N3+1];
 868                 double y4 = coord[dim_space*N4+1];
 869                 double y5 = coord[dim_space*N5+1];
 870                 double y6 = coord[dim_space*N6+1];
 871
 872                 double z1 = coord[dim_space*N1+2];
 873                 double z2 = coord[dim_space*N2+2];
 874                 double z3 = coord[dim_space*N3+2];
 875                 double z4 = coord[dim_space*N4+2];
 876                 double z5 = coord[dim_space*N5+2];
 877                 double z6 = coord[dim_space*N6+2];
 878
 879                 double a1 = (x2-x3)/2.0, a2 = (y2-y3)/2.0, a3 = (z2-z3)/2.0;
 880                 double b1 = (x5-x6)/2.0, b2 = (y5-y6)/2.0, b3 = (z5-z6)/2.0;
 881                 double c1 = (x4-x1)/2.0, c2 = (y4-y1)/2.0, c3 = (z4-z1)/2.0;
 882                 double d1 = (x5-x2)/2.0, d2 = (y5-y2)/2.0, d3 = (z5-z2)/2.0;
 883                 double e1 = (x6-x3)/2.0, e2 = (y6-y3)/2.0, e3 = (z6-z3)/2.0;
 884                 double f1 = (x1-x3)/2.0, f2 = (y1-y3)/2.0, f3 = (z1-z3)/2.0;
 885                 double h1 = (x4-x6)/2.0, h2 = (y4-y6)/2.0, h3 = (z4-z6)/2.0;
 886
 887                 double A = a1*c2*f3 - a1*c3*f2 - a2*c1*f3 + a2*c3*f1 +
 888                   a3*c1*f2 - a3*c2*f1;
 889                 double B = b1*c2*h3 - b1*c3*h2 - b2*c1*h3 + b2*c3*h1 +
 890                   b3*c1*h2 - b3*c2*h1;
 891                 double C = (a1*c2*h3 + b1*c2*f3) - (a1*c3*h2 + b1*c3*f2) -
 892                   (a2*c1*h3 + b2*c1*f3) + (a2*c3*h1 + b2*c3*f1) +
 893                   (a3*c1*h2 + b3*c1*f2) - (a3*c2*h1 + b3*c2*f1);
 894                 double D = a1*d2*f3 - a1*d3*f2 - a2*d1*f3 + a2*d3*f1 +
 895                   a3*d1*f2 - a3*d2*f1;
 896                 double E = b1*d2*h3 - b1*d3*h2 - b2*d1*h3 + b2*d3*h1 +
 897                   b3*d1*h2 - b3*d2*h1;
 898                 double F = (a1*d2*h3 + b1*d2*f3) - (a1*d3*h2 + b1*d3*f2) -
 899                   (a2*d1*h3 + b2*d1*f3) + (a2*d3*h1 + b2*d3*f1) +
 900                   (a3*d1*h2 + b3*d1*f2) - (a3*d2*h1 + b3*d2*f1);
 901                 double G = a1*e2*f3 - a1*e3*f2 - a2*e1*f3 + a2*e3*f1 +
 902                   a3*e1*f2 - a3*e2*f1;
 903                 double H = b1*e2*h3 - b1*e3*h2 - b2*e1*h3 + b2*e3*h1 +
 904                   b3*e1*h2 - b3*e2*h1;
 905                 double P = (a1*e2*h3 + b1*e2*f3) - (a1*e3*h2 + b1*e3*f2) -
 906                   (a2*e1*h3 + b2*e1*f3) + (a2*e3*h1 + b2*e3*f1) +
 907                   (a3*e1*h2 + b3*e1*f2) - (a3*e2*h1 + b3*e2*f1);
 908
 909                 xvolume = -2.0*(2.0*(A + B + D + E + G + H) + C + F + P)/9.0;
 910
 911                 //volume[index] = xvolume ;
 912                 volume->setIJ(index,1,xvolume) ;
 913                 index++;
 914               }
 915             break;
 916           }
 917         case MED_HEXA8 : case MED_HEXA20 :
 918           {
 919             for (int hexa=0;hexa<nb_entity_type;hexa++)
 920               {
 921                 int hexa_index = (type%100)*hexa;
 922
 923                 int N1 = global_connectivity[hexa_index]-1;
 924                 int N2 = global_connectivity[hexa_index+1]-1;
 925                 int N3 = global_connectivity[hexa_index+2]-1;
 926                 int N4 = global_connectivity[hexa_index+3]-1;
 927                 int N5 = global_connectivity[hexa_index+4]-1;
 928                 int N6 = global_connectivity[hexa_index+5]-1;
 929                 int N7 = global_connectivity[hexa_index+6]-1;
 930                 int N8 = global_connectivity[hexa_index+7]-1;
 931
 932                 double x1 = coord[dim_space*N1];
 933                 double x2 = coord[dim_space*N2];
 934                 double x3 = coord[dim_space*N3];
 935                 double x4 = coord[dim_space*N4];
 936                 double x5 = coord[dim_space*N5];
 937                 double x6 = coord[dim_space*N6];
 938                 double x7 = coord[dim_space*N7];
 939                 double x8 = coord[dim_space*N8];
 940
 941                 double y1 = coord[dim_space*N1+1];
 942                 double y2 = coord[dim_space*N2+1];
 943                 double y3 = coord[dim_space*N3+1];
 944                 double y4 = coord[dim_space*N4+1];
 945                 double y5 = coord[dim_space*N5+1];
 946                 double y6 = coord[dim_space*N6+1];
 947                 double y7 = coord[dim_space*N7+1];
 948                 double y8 = coord[dim_space*N8+1];
 949
 950                 double z1 = coord[dim_space*N1+2];
 951                 double z2 = coord[dim_space*N2+2];
 952                 double z3 = coord[dim_space*N3+2];
 953                 double z4 = coord[dim_space*N4+2];
 954                 double z5 = coord[dim_space*N5+2];
 955                 double z6 = coord[dim_space*N6+2];
 956                 double z7 = coord[dim_space*N7+2];
 957                 double z8 = coord[dim_space*N8+2];
 958
 959                 double a1 = (x3-x4)/8.0, a2 = (y3-y4)/8.0, a3 = (z3-z4)/8.0;
 960                 double b1 = (x2-x1)/8.0, b2 = (y2-y1)/8.0, b3 = (z2-z1)/8.0;
 961                 double c1 = (x7-x8)/8.0, c2 = (y7-y8)/8.0, c3 = (z7-z8)/8.0;
 962                 double d1 = (x6-x5)/8.0, d2 = (y6-y5)/8.0, d3 = (z6-z5)/8.0;
 963                 double e1 = (x3-x2)/8.0, e2 = (y3-y2)/8.0, e3 = (z3-z2)/8.0;
 964                 double f1 = (x4-x1)/8.0, f2 = (y4-y1)/8.0, f3 = (z4-z1)/8.0;
 965                 double h1 = (x7-x6)/8.0, h2 = (y7-y6)/8.0, h3 = (z7-z6)/8.0;
 966                 double p1 = (x8-x5)/8.0, p2 = (y8-y5)/8.0, p3 = (z8-z5)/8.0;
 967                 double q1 = (x3-x7)/8.0, q2 = (y3-y7)/8.0, q3 = (z3-z7)/8.0;
 968                 double r1 = (x4-x8)/8.0, r2 = (y4-y8)/8.0, r3 = (z4-z8)/8.0;
 969                 double s1 = (x2-x6)/8.0, s2 = (y2-y6)/8.0, s3 = (z2-z6)/8.0;
 970                 double t1 = (x1-x5)/8.0, t2 = (y1-y5)/8.0, t3 = (z1-z5)/8.0;
 971
 972                 double A = a1*e2*q3 - a1*e3*q2 - a2*e1*q3 + a2*e3*q1 +
 973                   a3*e1*q2 - a3*e2*q1;
 974                 double B = c1*h2*q3 - c1*h3*q2 - c2*h1*q3 + c2*h3*q1 +
 975                   c3*h1*q2 - c3*h2*q1;
 976                 double C = (a1*h2 + c1*e2)*q3 - (a1*h3 + c1*e3)*q2 -
 977                   (a2*h1 + c2*e1)*q3 + (a2*h3 + c2*e3)*q1 +
 978                   (a3*h1 + c3*e1)*q2 - (a3*h2 + c3*e2)*q1;
 979                 double D = b1*e2*s3 - b1*e3*s2 - b2*e1*s3 + b2*e3*s1 +
 980                   b3*e1*s2 - b3*e2*s1;
 981                 double E = d1*h2*s3 - d1*h3*s2 - d2*h1*s3 + d2*h3*s1 +
 982                   d3*h1*s2 - d3*h2*s1;
 983                 double F = (b1*h2 + d1*e2)*s3 - (b1*h3 + d1*e3)*s2 -
 984                   (b2*h1 + d2*e1)*s3 + (b2*h3 + d2*e3)*s1 +
 985                   (b3*h1 + d3*e1)*s2 - (b3*h2 + d3*e2)*s1;
 986                 double G = (a1*e2*s3 + b1*e2*q3) - (a1*e3*s2 + b1*e3*q2) -
 987                   (a2*e1*s3 + b2*e1*q3) + (a2*e3*s1 + b2*e3*q1) +
 988                   (a3*e1*s2 + b3*e1*q2) - (a3*e2*s1 + b3*e2*q1);
 989                 double H = (c1*h2*s3 + d1*h2*q3) - (c1*h3*s2 + d1*h3*q2) -
 990                   (c2*h1*s3 + d2*h1*q3) + (c2*h3*s1 + d2*h3*q1) +
 991                   (c3*h1*s2 + d3*h1*q2) - (c3*h2*s1 + d3*h2*q1);
 992                 double I = ((a1*h2 + c1*e2)*s3 + (b1*h2 + d1*e2)*q3) -
 993                   ((a1*h3 + c1*e3)*s2 + (b1*h3 + d1*e3)*q2) -
 994                   ((a2*h1 + c2*e1)*s3 + (b2*h1 + d2*e1)*q3) +
 995                   ((a2*h3 + c2*e3)*s1 + (b2*h3 + d2*e3)*q1) +
 996                   ((a3*h1 + c3*e1)*s2 + (b3*h1 + d3*e1)*q2) -
 997                   ((a3*h2 + c3*e2)*s1 + (b3*h2 + d3*e2)*q1);
 998                 double J = a1*f2*r3 - a1*f3*r2 - a2*f1*r3 + a2*f3*r1 +
 999                   a3*f1*r2 - a3*f2*r1;
1000                 double K = c1*p2*r3 - c1*p3*r2 - c2*p1*r3 + c2*p3*r1 +
1001                   c3*p1*r2 - c3*p2*r1;
1002                 double L = (a1*p2 + c1*f2)*r3 - (a1*p3 + c1*f3)*r2 -
1003                   (a2*p1 + c2*f1)*r3 + (a2*p3 + c2*f3)*r1 +
1004                   (a3*p1 + c3*f1)*r2 - (a3*p2 + c3*f2)*r1;
1005                 double M = b1*f2*t3 - b1*f3*t2 - b2*f1*t3 + b2*f3*t1 +
1006                   b3*f1*t2 - b3*f2*t1;
1007                 double N = d1*p2*t3 - d1*p3*t2 - d2*p1*t3 + d2*p3*t1 +
1008                   d3*p1*t2 - d3*p2*t1;
1009                 double O = (b1*p2 + d1*f2)*t3 - (b1*p3 + d1*f3)*t2 -
1010                   (b2*p1 + d2*f1)*t3 + (b2*p3 + d2*f3)*t1 +
1011                   (b3*p1 + d3*f1)*t2 - (b3*p2 + d3*f2)*t1;
1012                 double P = (a1*f2*t3 + b1*f2*r3) - (a1*f3*t2 + b1*f3*r2) -
1013                   (a2*f1*t3 + b2*f1*r3) + (a2*f3*t1 + b2*f3*r1) +
1014                   (a3*f1*t2 + b3*f1*r2) - (a3*f2*t1 + b3*f2*r1);
1015                 double Q = (c1*p2*t3 + d1*p2*r3) - (c1*p3*t2 + d1*p3*r2) -
1016                   (c2*p1*t3 + d2*p1*r3) + (c2*p3*t1 + d2*p3*r1) +
1017                   (c3*p1*t2 + d3*p1*r2) - (c3*p2*t1 + d3*p2*r1);
1018                 double R = ((a1*p2 + c1*f2)*t3 + (b1*p2 + d1*f2)*r3) -
1019                   ((a1*p3 + c1*f3)*t2 + (b1*p3 + d1*f3)*r2) -
1020                   ((a2*p1 + c2*f1)*t3 + (b2*p1 + d2*f1)*r3) +
1021                   ((a2*p3 + c2*f3)*t1 + (b2*p3 + d2*f3)*r1) +
1022                   ((a3*p1 + c3*f1)*t2 + (b3*p1 + d3*f1)*r2) -
1023                   ((a3*p2 + c3*f2)*t1 + (b3*p2 + d3*f2)*r1);
1024                 double S = (a1*e2*r3 + a1*f2*q3) - (a1*e3*r2 + a1*f3*q2) -
1025                   (a2*e1*r3 + a2*f1*q3) + (a2*e3*r1 + a2*f3*q1) +
1026                   (a3*e1*r2 + a3*f1*q2) - (a3*e2*r1 + a3*f2*q1);
1027                 double T = (c1*h2*r3 + c1*p2*q3) - (c1*h3*r2 + c1*p3*q2) -
1028                   (c2*h1*r3 + c2*p1*q3) + (c2*h3*r1 + c2*p3*q1) +
1029                   (c3*h1*r2 + c3*p1*q2) - (c3*h2*r1 + c3*p2*q1);
1030                 double U = ((a1*h2 + c1*e2)*r3 + (a1*p2 + c1*f2)*q3) -
1031                   ((a1*h3 + c1*e3)*r2 + (a1*p3 + c1*f3)*q2) -
1032                   ((a2*h1 + c2*e1)*r3 + (a2*p1 + c2*f1)*q3) +
1033                   ((a2*h3 + c2*e3)*r1 + (a2*p3 + c2*f3)*q1) +
1034                   ((a3*h1 + c3*e1)*r2 + (a3*p1 + c3*f1)*q2) -
1035                   ((a3*h2 + c3*e2)*r1 + (a3*p2 + c3*f2)*q1);
1036                 double V = (b1*e2*t3 + b1*f2*s3) - (b1*e3*t2 + b1*f3*s2) -
1037                   (b2*e1*t3 + b2*f1*s3) + (b2*e3*t1 + b2*f3*s1) +
1038                   (b3*e1*t2 + b3*f1*s2) - (b3*e2*t1 + b3*f2*s1);
1039                 double W = (d1*h2*t3 + d1*p2*s3) - (d1*h3*t2 + d1*p3*s2) -
1040                   (d2*h1*t3 + d2*p1*s3) + (d2*h3*t1 + d2*p3*s1) +
1041                   (d3*h1*t2 + d3*p1*s2) - (d3*h2*t1 + d3*p2*s1);
1042                 double X = ((b1*h2 + d1*e2)*t3 + (b1*p2 + d1*f2)*s3) -
1043                   ((b1*h3 + d1*e3)*t2 + (b1*p3 + d1*f3)*s2) -
1044                   ((b2*h1 + d2*e1)*t3 + (b2*p1 + d2*f1)*s3) +
1045                   ((b2*h3 + d2*e3)*t1 + (b2*p3 + d2*f3)*s1) +
1046                   ((b3*h1 + d3*e1)*t2 + (b3*p1 + d3*f1)*s2) -
1047                   ((b3*h2 + d3*e2)*t1 + (b3*p2 + d3*f2)*s1);
1048                 double Y = (a1*e2*t3 + a1*f2*s3 + b1*e2*r3 + b1*f2*q3) -
1049                   (a1*e3*t2 + a1*f3*s2 + b1*e3*r2 + b1*f3*q2) -
1050                   (a2*e1*t3 + a2*f1*s3 + b2*e1*r3 + b2*f1*q3) +
1051                   (a2*e3*t1 + a2*f3*s1 + b2*e3*r1 + b2*f3*q1) +
1052                   (a3*e1*t2 + a3*f1*s2 + b3*e1*r2 + b3*f1*q2) -
1053                   (a3*e2*t1 + a3*f2*s1 + b3*e2*r1 + b3*f2*q1);
1054                 double Z = (c1*h2*t3 + c1*p2*s3 + d1*h2*r3 + d1*p2*q3) -
1055                   (c1*h3*t2 + c1*p3*s2 + d1*h3*r2 + d1*p3*q2) -
1056                   (c2*h1*t3 + c2*p1*s3 + d2*h1*r3 + d2*p1*q3) +
1057                   (c2*h3*t1 + c2*p3*s1 + d2*h3*r1 + d2*p3*q1) +
1058                   (c3*h1*t2 + c3*p1*s2 + d3*h1*r2 + d3*p1*q2) -
1059                   (c3*h2*t1 + c3*p2*s1 + d3*h2*r1 + d3*p2*q1);
1060                 double AA = ((a1*h2 + c1*e2)*t3 + (a1*p2 + c1*f2)*s3 +
1061                              (b1*h2 + d1*e2)*r3 + (b1*p2 + d1*f2)*q3) -
1062                   ((a1*h3 + c1*e3)*t2 + (a1*p3 + c1*f3)*s2 +
1063                    (b1*h3 + d1*e3)*r2 + (b1*p3 + d1*f3)*q2) -
1064                   ((a2*h1 + c2*e1)*t3 + (a2*p1 + c2*f1)*s3 +
1065                    (b2*h1 + d2*e1)*r3 + (b2*p1 + d2*f1)*q3) +
1066                   ((a2*h3 + c2*e3)*t1 + (a2*p3 + c2*f3)*s1 +
1067                    (b2*h3 + d2*e3)*r1 + (b2*p3 + d2*f3)*q1) +
1068                   ((a3*h1 + c3*e1)*t2 + (a3*p1 + c3*f1)*s2 +
1069                    (b3*h1 + d3*e1)*r2 + (b3*p1 + d3*f1)*q2) -
1070                   ((a3*h2 + c3*e2)*t1 + (a3*p2 + c3*f2)*s1 +
1071                    (b3*h2 + d3*e2)*r1 + (b3*p2 + d3*f2)*q1);
1072
1073                 xvolume = 64.0*(8.0*(A + B + D + E + J + K + M + N) +
1074                                 4.0*(C + F + G + H + L + O + P + Q + S + T +
1075                                      V + W) + 2.0*(I + R + U + X + Y + Z) +
1076                                 AA)/27.0;
1077
1078                 //volume[index] = xvolume ;
1079                 volume->setIJ(index,1,xvolume) ;
1080                 index++;
1081               }
1082             break;
1083           }
1084         default :
1085           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Volumes on it !"));
1086           break;
1087         }
1088     }
1089
1090   return Volume;
1091 }
1092
1093 FIELD<double>* MESH::getArea(const SUPPORT * Support) const throw (MEDEXCEPTION)
1094 {
1095   const char * LOC = "MESH::getArea(SUPPORT*) : ";
1096   BEGIN_OF(LOC);
1097
1098   // Support must be on 2D elements
1099
1100   // Make sure that the MESH class is the same as the MESH class attribut
1101   // in the class Support
1102   MESH* myMesh = Support->getMesh();
1103   if (this != myMesh)
1104     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1105
1106   int dim_space = getSpaceDimension();
1107   medEntityMesh support_entity = Support->getEntity();
1108   bool onAll = Support->isOnAllElements();
1109
1110   int nb_type, length_values;
1111   const medGeometryElement* types;
1112   int nb_entity_type;
1113   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1114   const int* global_connectivity;
1115
1116 //    if (onAll)
1117 //      {
1118 //        nb_type = myMesh->getNumberOfTypes(support_entity);
1119 //        length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1120 //        types = getTypes(support_entity);
1121 //      }
1122 //    else
1123     {
1124       nb_type = Support->getNumberOfTypes();
1125       length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1126       types = Support->getTypes();
1127     }
1128
1129   int index;
1130   FIELD<double>* Area;
1131
1132   Area = new FIELD<double>::FIELD(Support,1);
1133   Area->setName("AREA");
1134   Area->setDescription("cells or faces area");
1135
1136   Area->setComponentName(1,"area");
1137   Area->setComponentDescription(1,"desc-comp");
1138
1139   /*  string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1];*/
1140
1141   string MEDComponentUnit="trucmuch";
1142
1143   Area->setMEDComponentUnit(1,MEDComponentUnit);
1144
1145   Area->setValueType(MED_REEL64);
1146
1147   Area->setIterationNumber(0);
1148   Area->setOrderNumber(0);
1149   Area->setTime(0.0);
1150
1151   double *area = new double[length_values];
1152   //double *area = Area->getValue(MED_FULL_INTERLACE);
1153
1154   const double * coord = getCoordinates(MED_FULL_INTERLACE);
1155   index = 0;
1156
1157   for (int i=0;i<nb_type;i++)
1158     {
1159       medGeometryElement type = types[i] ;
1160       double xarea;
1161
1162       if (onAll)
1163         {
1164           nb_entity_type = getNumberOfElements(support_entity,type);
1165           global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1166         }
1167       else
1168         {
1169           // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all !"));
1170
1171           nb_entity_type = Support->getNumberOfElements(type);
1172
1173           const int * supp_number = Support->getNumber(type);
1174           const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1175           const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1176           int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1177
1178           for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1179             for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1180               global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1181             }
1182           }
1183
1184           global_connectivity = global_connectivity_tmp ;
1185
1186         }
1187
1188       switch (type)
1189         {
1190         case MED_TRIA3 : case MED_TRIA6 :
1191           {
1192             for (int tria=0;tria<nb_entity_type;tria++)
1193               {
1194                 int tria_index = (type%100)*tria;
1195
1196                 int N1 = global_connectivity[tria_index]-1;
1197                 int N2 = global_connectivity[tria_index+1]-1;
1198                 int N3 = global_connectivity[tria_index+2]-1;
1199
1200                 double x1 = coord[dim_space*N1];
1201                 double x2 = coord[dim_space*N2];
1202                 double x3 = coord[dim_space*N3];
1203
1204                 double y1 = coord[dim_space*N1+1];
1205                 double y2 = coord[dim_space*N2+1];
1206                 double y3 = coord[dim_space*N3+1];
1207
1208                 if (dim_space==2)
1209                   {
1210                     xarea = - ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))/2.0;
1211                   }
1212                 else
1213                   {
1214                     double z1 = coord[dim_space*N1+2];
1215                     double z2 = coord[dim_space*N2+2];
1216                     double z3 = coord[dim_space*N3+2];
1217
1218                     xarea = sqrt(((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1))*
1219                                  ((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1)) +
1220                                  ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1))*
1221                                  ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1)) +
1222                                  ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))*
1223                                  ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1)))/2.0;
1224                   }
1225
1226                 area[index] = xarea ;
1227                 index++;
1228               }
1229             break;
1230           }
1231         case MED_QUAD4 : case MED_QUAD8 :
1232           {
1233             for (int quad=0;quad<nb_entity_type;quad++)
1234               {
1235                 int quad_index = (type%100)*quad;
1236
1237                 int N1 = global_connectivity[quad_index]-1;
1238                 int N2 = global_connectivity[quad_index+1]-1;
1239                 int N3 = global_connectivity[quad_index+2]-1;
1240                 int N4 = global_connectivity[quad_index+3]-1;
1241
1242                 double x1 = coord[dim_space*N1];
1243                 double x2 = coord[dim_space*N2];
1244                 double x3 = coord[dim_space*N3];
1245                 double x4 = coord[dim_space*N4];
1246
1247                 double y1 = coord[dim_space*N1+1];
1248                 double y2 = coord[dim_space*N2+1];
1249                 double y3 = coord[dim_space*N3+1];
1250                 double y4 = coord[dim_space*N4+1];
1251
1252                 if (dim_space==2)
1253                   {
1254                     double a1 = (x2-x1)/4.0, a2 = (y2-y1)/4.0;
1255                     double b1 = (x3-x4)/4.0, b2 = (y3-y4)/4.0;
1256                     double c1 = (x3-x2)/4.0, c2 = (y3-y2)/4.0;
1257                     double d1 = (x4-x1)/4.0, d2 = (y4-y1)/4.0;
1258
1259                     xarea = - 4.0*(b1*c2 - c1*b2 + a1*c2 - c1*a2 + b1*d2 -
1260                                  d1*b2 + a1*d2 - d1*a2);
1261                   }
1262                 else
1263                   {
1264                     double z1 = coord[dim_space*N1+2];
1265                     double z2 = coord[dim_space*N2+2];
1266                     double z3 = coord[dim_space*N3+2];
1267                     double z4 = coord[dim_space*N4+2];
1268
1269                     xarea = (sqrt(((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1))*
1270                                   ((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1)) +
1271                                   ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1))*
1272                                   ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1)) +
1273                                   ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))*
1274                                   ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))) +
1275                              sqrt(((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3))*
1276                                   ((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3)) +
1277                                   ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3))*
1278                                   ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3)) +
1279                                   ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))*
1280                                   ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))))/2.0;
1281                   }
1282
1283                 area[index] = xarea ;
1284                 index++;
1285               }
1286             break;
1287           }
1288         default :
1289           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Areas on it !"));
1290           break;
1291         }
1292     }
1293
1294   Area->setValue(MED_FULL_INTERLACE,area);
1295   delete[] area;
1296   return Area;
1297 }
1298
1299 FIELD<double>* MESH::getLength(const SUPPORT * Support) const throw (MEDEXCEPTION)
1300 {
1301   const char * LOC = "MESH::getLength(SUPPORT*) : ";
1302   BEGIN_OF(LOC);
1303
1304   // Support must be on 1D elements
1305
1306   // Make sure that the MESH class is the same as the MESH class attribut
1307   // in the class Support
1308   MESH* myMesh = Support->getMesh();
1309   if (this != myMesh)
1310     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1311
1312   int dim_space = getSpaceDimension();
1313   medEntityMesh support_entity = Support->getEntity();
1314   bool onAll = Support->isOnAllElements();
1315
1316   int nb_type, length_values;
1317   const medGeometryElement* types;
1318   int nb_entity_type;
1319   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1320   const int* global_connectivity;
1321
1322 //    if (onAll)
1323 //      {
1324 //        nb_type = myMesh->getNumberOfTypes(support_entity);
1325 //        length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1326 //        types = getTypes(support_entity);
1327 //      }
1328 //    else
1329     {
1330       nb_type = Support->getNumberOfTypes();
1331       length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1332       types = Support->getTypes();
1333     }
1334
1335   int index;
1336   FIELD<double>* Length;
1337
1338   Length = new FIELD<double>::FIELD(Support,1);
1339   //  double *length = new double [length_values];
1340   Length->setValueType(MED_REEL64);
1341
1342   //const double *length = Length->getValue(MED_FULL_INTERLACE);
1343   MEDARRAY<double> * length = Length->getvalue();
1344
1345   const double * coord = getCoordinates(MED_FULL_INTERLACE);
1346   index = 1;
1347
1348   for (int i=0;i<nb_type;i++)
1349     {
1350       medGeometryElement type = types[i] ;
1351       double xlength;
1352
1353       if (onAll)
1354         {
1355           nb_entity_type = getNumberOfElements(support_entity,type);
1356           global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1357         }
1358       else
1359         {
1360           //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all !"));
1361
1362           nb_entity_type = Support->getNumberOfElements(type);
1363
1364           const int * supp_number = Support->getNumber(type);
1365           const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1366           const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1367           int* global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1368
1369           for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1370             for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1371               global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1372             }
1373           }
1374           global_connectivity = global_connectivity_tmp ;
1375         }
1376
1377       switch (type)
1378         {
1379         case MED_SEG2 : case MED_SEG3 :
1380           {
1381             for (int edge=0;edge<nb_entity_type;edge++)
1382               {
1383                 int edge_index = (type%100)*edge;
1384
1385                 int N1 = global_connectivity[edge_index]-1;
1386                 int N2 = global_connectivity[edge_index+1]-1;
1387
1388                 double x1 = coord[dim_space*N1];
1389                 double x2 = coord[dim_space*N2];
1390
1391                 double y1 = coord[dim_space*N1+1];
1392                 double y2 = coord[dim_space*N2+1];
1393
1394                 double z1, z2 ; z1 = 0.0 ; z2 = 0.0 ; if (dim_space==3) {
1395                   z1 = coord[dim_space*N1+2]; z2 = coord[dim_space*N2+2];}
1396
1397                 xlength =  sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) +
1398                                 (z1 - z2)*(z1 - z2));
1399
1400                 length->setIJ(index,1,xlength) ;
1401                 index++;
1402               }
1403             break;
1404           }
1405         default :
1406           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Lengths on it !"));
1407           break;
1408         }
1409     }
1410
1411   return Length;
1412 }
1413
1414 FIELD<double>* MESH::getNormal(const SUPPORT * Support) const throw (MEDEXCEPTION)
1415 {
1416   const char * LOC = "MESH::getNormal(SUPPORT*) : ";
1417   BEGIN_OF(LOC);
1418
1419   // Support must be on 2D or 1D elements
1420
1421   // Make sure that the MESH class is the same as the MESH class attribut
1422   // in the class Support
1423   MESH* myMesh = Support->getMesh();
1424   if (this != myMesh)
1425     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : pointeur problem !"));
1426
1427   int dim_space = getSpaceDimension();
1428   medEntityMesh support_entity = Support->getEntity();
1429   bool onAll = Support->isOnAllElements();
1430
1431   int nb_type, length_values;
1432   const medGeometryElement* types;
1433   int nb_entity_type;
1434   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1435   const int* global_connectivity;
1436
1437 //    if (onAll)
1438 //      {
1439 //        nb_type = myMesh->getNumberOfTypes(support_entity);
1440 //        length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1441 //        types = getTypes(support_entity);
1442 //      }
1443 //    else
1444     {
1445       nb_type = Support->getNumberOfTypes();
1446       length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1447       types = Support->getTypes();
1448     }
1449
1450   int index;
1451
1452   FIELD<double>* Normal = new FIELD<double>::FIELD(Support,dim_space);
1453   Normal->setName("NORMAL");
1454   Normal->setDescription("cells or faces normal");
1455   for (int k=1;k<=dim_space;k++) {
1456     Normal->setComponentName(k,"normal");
1457     Normal->setComponentDescription(k,"desc-comp");
1458     Normal->setMEDComponentUnit(k,"unit");
1459   }
1460
1461   Normal->setValueType(MED_REEL64);
1462
1463   Normal->setIterationNumber(MED_NOPDT);
1464   Normal->setOrderNumber(MED_NONOR);
1465   Normal->setTime(0.0);
1466
1467   double * normal = new double [dim_space*length_values];
1468   //double *normal = Normal->getValue(MED_FULL_INTERLACE);
1469
1470   const double * coord = getCoordinates(MED_FULL_INTERLACE);
1471   index = 0;
1472
1473   for (int i=0;i<nb_type;i++)
1474     {
1475       medGeometryElement type = types[i] ;
1476
1477       // Make sure we trying to get Normals on
1478       // TRIA3 or TRIA6 or QUAD4 QUAD8 (2D) cells on a 3D mesh
1479       // or on SEG2 or SEG3 (1D) cells on a 2D mesh
1480
1481       if ( (((type==MED_TRIA3) || (type==MED_TRIA6) ||
1482              (type==MED_QUAD4) || (type==MED_QUAD8)) &&
1483             (dim_space != 3)) || (((type==MED_SEG2) || (type==MED_SEG3)) &&
1484                                   (dim_space != 2)) )
1485         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : dimension problem !"));
1486
1487       double xnormal1, xnormal2, xnormal3 ;
1488
1489       if (onAll)
1490         {
1491           nb_entity_type = getNumberOfElements(support_entity,type);
1492           global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1493         }
1494       else
1495         {
1496           // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all for instance !"));
1497           nb_entity_type = Support->getNumberOfElements(type);
1498
1499           const int * supp_number = Support->getNumber(type);
1500           const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1501           const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1502           int * global_connectivity_tmp = 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_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1507             }
1508           }
1509
1510           global_connectivity = global_connectivity_tmp ;
1511
1512         }
1513
1514       switch (type)
1515         {
1516         case MED_TRIA3 : case MED_TRIA6 :
1517           {
1518             for (int tria=0;tria<nb_entity_type;tria++)
1519               {
1520                 int tria_index = (type%100)*tria;
1521
1522                 int N1 = global_connectivity[tria_index]-1;
1523                 int N2 = global_connectivity[tria_index+1]-1;
1524                 int N3 = global_connectivity[tria_index+2]-1;
1525
1526                 double xarea;
1527                 double x1 = coord[dim_space*N1];
1528                 double x2 = coord[dim_space*N2];
1529                 double x3 = coord[dim_space*N3];
1530
1531                 double y1 = coord[dim_space*N1+1];
1532                 double y2 = coord[dim_space*N2+1];
1533                 double y3 = coord[dim_space*N3+1];
1534
1535                 double z1 = coord[dim_space*N1+2];
1536                 double z2 = coord[dim_space*N2+2];
1537                 double z3 = coord[dim_space*N3+2];
1538
1539                 xnormal1 = ((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1))/2.0;
1540                 xnormal2 = ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1))/2.0;
1541                 xnormal3 = ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))/2.0;
1542
1543                 normal[3*index] = xnormal1 ;
1544                 normal[3*index+1] = xnormal2 ;
1545                 normal[3*index+2] = xnormal3 ;
1546
1547                 index++;
1548               }
1549             break;
1550           }
1551         case MED_QUAD4 : case MED_QUAD8 :
1552           {
1553             for (int quad=0;quad<nb_entity_type;quad++)
1554               {
1555                 int quad_index = (type%100)*quad;
1556
1557                 int N1 = global_connectivity[quad_index]-1;
1558                 int N2 = global_connectivity[quad_index+1]-1;
1559                 int N3 = global_connectivity[quad_index+2]-1;
1560                 int N4 = global_connectivity[quad_index+3]-1;
1561
1562                 double xarea;
1563                 double x1 = coord[dim_space*N1];
1564                 double x2 = coord[dim_space*N2];
1565                 double x3 = coord[dim_space*N3];
1566                 double x4 = coord[dim_space*N4];
1567
1568                 double y1 = coord[dim_space*N1+1];
1569                 double y2 = coord[dim_space*N2+1];
1570                 double y3 = coord[dim_space*N3+1];
1571                 double y4 = coord[dim_space*N4+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
1578                 xnormal1 = (y2-y1)*(z4-z1) - (y4-y1)*(z2-z1);
1579                 xnormal2 = (x4-x1)*(z2-z1) - (x2-x1)*(z4-z1);
1580                 xnormal3 = (x2-x1)*(y4-y1) - (x4-x1)*(y2-y1);
1581                 xarea = sqrt(xnormal1*xnormal1 + xnormal2*xnormal2 +
1582                              xnormal3*xnormal3);
1583
1584                 xnormal1 = xnormal1/xarea;
1585                 xnormal2 = xnormal2/xarea;
1586                 xnormal3 = xnormal3/xarea;
1587
1588                 xarea = (sqrt(((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1))*
1589                               ((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1)) +
1590                               ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1))*
1591                               ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1)) +
1592                               ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))*
1593                               ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))) +
1594                          sqrt(((y3-y2)*(z4-z3) - (y4-y3)*(z3-z2))*
1595                               ((y3-y2)*(z4-z3) - (y4-y3)*(z3-z2)) +
1596                               ((x4-x3)*(z2-z2) - (x3-x2)*(z4-z3))*
1597                               ((x4-x3)*(z4-z2) - (x3-x2)*(z4-z3)) +
1598                               ((x3-x2)*(y4-y3) - (x4-x3)*(y3-y2))*
1599                               ((x3-x2)*(y4-y3) - (x4-x3)*(y3-y2))))/2.0;
1600
1601                 xnormal1 = xnormal1*xarea;
1602                 xnormal2 = xnormal2*xarea;
1603                 xnormal3 = xnormal3*xarea;
1604
1605                 normal[3*index] = xnormal1 ;
1606                 normal[3*index+1] = xnormal2 ;
1607                 normal[3*index+2] = xnormal3 ;
1608
1609                 index++;
1610               }
1611             break;
1612           }
1613         case MED_SEG2 : case MED_SEG3 :
1614           {
1615             for (int edge=0;edge<nb_entity_type;edge++)
1616               {
1617                 int edge_index = (type%100)*edge;
1618
1619                 int N1 = global_connectivity[edge_index]-1;
1620                 int N2 = global_connectivity[edge_index+1]-1;
1621
1622                 double x1 = coord[dim_space*N1];
1623                 double x2 = coord[dim_space*N2];
1624
1625                 double y1 = coord[dim_space*N1+1];
1626                 double y2 = coord[dim_space*N2+1];
1627
1628                 xnormal1 = -(y2-y1);
1629                 xnormal2 = x2-x1;
1630
1631                 normal[2*index] = xnormal1 ;
1632                 normal[2*index+1] = xnormal2 ;
1633
1634                 index++;
1635               }
1636             break;
1637           }
1638         default :
1639           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Normals on it !"));
1640           break;
1641         }
1642     }
1643
1644   Normal->setValue(MED_FULL_INTERLACE,normal);
1645   delete[] normal ;
1646
1647   END_OF(LOC);
1648
1649   return Normal;
1650 }
1651
1652 FIELD<double>* MESH::getBarycenter(const SUPPORT * Support) const throw (MEDEXCEPTION)
1653 {
1654   const char * LOC = "MESH::getBarycenter(SUPPORT*) : ";
1655   BEGIN_OF(LOC);
1656
1657   // Make sure that the MESH class is the same as the MESH class attribut
1658   // in the class Support
1659   MESH* myMesh = Support->getMesh();
1660   if (this != myMesh)
1661     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1662
1663   int dim_space = getSpaceDimension();
1664   medEntityMesh support_entity = Support->getEntity();
1665   bool onAll = Support->isOnAllElements();
1666
1667   int nb_type, length_values;
1668   const medGeometryElement* types;
1669   int nb_entity_type;
1670   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1671   const int* global_connectivity;
1672
1673 //    if (onAll)
1674 //      {
1675 //        nb_type = myMesh->getNumberOfTypes(support_entity);
1676 //        length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1677 //        types = getTypes(support_entity);
1678 //      }
1679 //    else
1680     {
1681       nb_type = Support->getNumberOfTypes();
1682       length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1683       types = Support->getTypes();
1684     }
1685
1686   int index;
1687   FIELD<double>* Barycenter;
1688
1689   Barycenter = new FIELD<double>::FIELD(Support,dim_space);
1690   Barycenter->setName("BARYCENTER");
1691   Barycenter->setDescription("cells or faces barycenter");
1692
1693   for (int k=0;k<dim_space;k++) {
1694     int kp1 = k + 1;
1695     Barycenter->setComponentName(kp1,myMesh->getCoordinatesNames()[k]);
1696     Barycenter->setComponentDescription(kp1,"desc-comp");
1697     Barycenter->setMEDComponentUnit(kp1,myMesh->getCoordinatesUnits()[k]);
1698   }
1699
1700   Barycenter->setValueType(MED_REEL64);
1701
1702   Barycenter->setIterationNumber(0);
1703   Barycenter->setOrderNumber(0);
1704   Barycenter->setTime(0.0);
1705
1706   double *barycenter = new double [dim_space*length_values];
1707   //  double *barycenter = Barycenter->getValue(MED_FULL_INTERLACE);
1708
1709   const double * coord = getCoordinates(MED_FULL_INTERLACE);
1710   index = 0;
1711
1712   for (int i=0;i<nb_type;i++)
1713     {
1714       medGeometryElement type = types[i] ;
1715       double xbarycenter1, xbarycenter2, xbarycenter3;
1716
1717       if (onAll)
1718         {
1719           nb_entity_type = getNumberOfElements(support_entity,type);
1720           global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1721         }
1722       else
1723         {
1724           // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all !"));
1725           nb_entity_type = Support->getNumberOfElements(type);
1726
1727           const int * supp_number = Support->getNumber(type);
1728           const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1729           const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1730           int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1731
1732           for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1733             for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1734               global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1735             }
1736           }
1737           global_connectivity = global_connectivity_tmp;
1738         }
1739
1740       switch (type)
1741         {
1742         case MED_TETRA4 : case MED_TETRA10 :
1743           {
1744             for (int tetra =0;tetra<nb_entity_type;tetra++)
1745               {
1746                 int tetra_index = (type%100)*tetra;
1747
1748                 int N1 = global_connectivity[tetra_index]-1;
1749                 int N2 = global_connectivity[tetra_index+1]-1;
1750                 int N3 = global_connectivity[tetra_index+2]-1;
1751                 int N4 = global_connectivity[tetra_index+3]-1;
1752
1753                 double x1 = coord[dim_space*N1];
1754                 double x2 = coord[dim_space*N2];
1755                 double x3 = coord[dim_space*N3];
1756                 double x4 = coord[dim_space*N4];
1757
1758                 double y1 = coord[dim_space*N1+1];
1759                 double y2 = coord[dim_space*N2+1];
1760                 double y3 = coord[dim_space*N3+1];
1761                 double y4 = coord[dim_space*N4+1];
1762
1763                 double z1 = coord[dim_space*N1+2];
1764                 double z2 = coord[dim_space*N2+2];
1765                 double z3 = coord[dim_space*N3+2];
1766                 double z4 = coord[dim_space*N4+2];
1767
1768                 xbarycenter1 = (x1 + x2 + x3 + x4)/4.0;
1769                 xbarycenter2 = (y1 + y2 + y3 + y4)/4.0;
1770                 xbarycenter3 = (z1 + z2 + z3 + z4)/4.0;
1771                 barycenter[3*index] = xbarycenter1 ;
1772                 barycenter[3*index+1] = xbarycenter2 ;
1773                 barycenter[3*index+2] = xbarycenter3 ;
1774                 index++;
1775               }
1776             break;
1777           }
1778         case MED_PYRA5 : case MED_PYRA13 :
1779           {
1780             for (int pyra=0;pyra<nb_entity_type;pyra++)
1781               {
1782                 int pyra_index = (type%100)*pyra;
1783
1784                 int N1 = global_connectivity[pyra_index]-1;
1785                 int N2 = global_connectivity[pyra_index+1]-1;
1786                 int N3 = global_connectivity[pyra_index+2]-1;
1787                 int N4 = global_connectivity[pyra_index+3]-1;
1788                 int N5 = global_connectivity[pyra_index+4]-1;
1789
1790                 double x1 = coord[dim_space*N1];
1791                 double x2 = coord[dim_space*N2];
1792                 double x3 = coord[dim_space*N3];
1793                 double x4 = coord[dim_space*N4];
1794                 double x5 = coord[dim_space*N5];
1795
1796                 double y1 = coord[dim_space*N1+1];
1797                 double y2 = coord[dim_space*N2+1];
1798                 double y3 = coord[dim_space*N3+1];
1799                 double y4 = coord[dim_space*N4+1];
1800                 double y5 = coord[dim_space*N5+1];
1801
1802                 double z1 = coord[dim_space*N1+2];
1803                 double z2 = coord[dim_space*N2+2];
1804                 double z3 = coord[dim_space*N3+2];
1805                 double z4 = coord[dim_space*N4+2];
1806                 double z5 = coord[dim_space*N5+2];
1807
1808                 xbarycenter1 = (x1 + x2 + x3 + x4 + x5)/5.0;
1809                 xbarycenter2 = (y1 + y2 + y3 + y4 + y5)/5.0;
1810                 xbarycenter3 = (z1 + z2 + z3 + z4 + z5)/5.0;
1811                 barycenter[3*index] = xbarycenter1 ;
1812                 barycenter[3*index+1] = xbarycenter2 ;
1813                 barycenter[3*index+2] = xbarycenter3 ;
1814                 index++;
1815               }
1816             break;
1817           }
1818         case MED_PENTA6 : case MED_PENTA15 :
1819           {
1820             for (int penta=0;penta<nb_entity_type;penta++)
1821               {
1822                 int penta_index = (type%100)*penta;
1823
1824                 int N1 = global_connectivity[penta_index]-1;
1825                 int N2 = global_connectivity[penta_index+1]-1;
1826                 int N3 = global_connectivity[penta_index+2]-1;
1827                 int N4 = global_connectivity[penta_index+3]-1;
1828                 int N5 = global_connectivity[penta_index+4]-1;
1829                 int N6 = global_connectivity[penta_index+5]-1;
1830
1831                 double x1 = coord[dim_space*N1];
1832                 double x2 = coord[dim_space*N2];
1833                 double x3 = coord[dim_space*N3];
1834                 double x4 = coord[dim_space*N4];
1835                 double x5 = coord[dim_space*N5];
1836                 double x6 = coord[dim_space*N6];
1837
1838                 double y1 = coord[dim_space*N1+1];
1839                 double y2 = coord[dim_space*N2+1];
1840                 double y3 = coord[dim_space*N3+1];
1841                 double y4 = coord[dim_space*N4+1];
1842                 double y5 = coord[dim_space*N5+1];
1843                 double y6 = coord[dim_space*N6+1];
1844
1845                 double z1 = coord[dim_space*N1+2];
1846                 double z2 = coord[dim_space*N2+2];
1847                 double z3 = coord[dim_space*N3+2];
1848                 double z4 = coord[dim_space*N4+2];
1849                 double z5 = coord[dim_space*N5+2];
1850                 double z6 = coord[dim_space*N6+2];
1851
1852                 xbarycenter1 = (x1 + x2 + x3 + x4 + x5 + x6)/6.0;
1853                 xbarycenter2 = (y1 + y2 + y3 + y4 + y5 + y6)/6.0;
1854                 xbarycenter3 = (z1 + z2 + z3 + z4 + z5 + z6)/6.0;
1855                 barycenter[3*index] = xbarycenter1 ;
1856                 barycenter[3*index+1] = xbarycenter2 ;
1857                 barycenter[3*index+2] = xbarycenter3 ;
1858                 index++;
1859               }
1860             break;
1861           }
1862         case MED_HEXA8 : case MED_HEXA20 :
1863           {
1864             for (int hexa=0;hexa<nb_entity_type;hexa++)
1865               {
1866                 int hexa_index = (type%100)*hexa;
1867
1868                 int N1 = global_connectivity[hexa_index]-1;
1869                 int N2 = global_connectivity[hexa_index+1]-1;
1870                 int N3 = global_connectivity[hexa_index+2]-1;
1871                 int N4 = global_connectivity[hexa_index+3]-1;
1872                 int N5 = global_connectivity[hexa_index+4]-1;
1873                 int N6 = global_connectivity[hexa_index+5]-1;
1874                 int N7 = global_connectivity[hexa_index+6]-1;
1875                 int N8 = global_connectivity[hexa_index+7]-1;
1876
1877                 double x1 = coord[dim_space*N1];
1878                 double x2 = coord[dim_space*N2];
1879                 double x3 = coord[dim_space*N3];
1880                 double x4 = coord[dim_space*N4];
1881                 double x5 = coord[dim_space*N5];
1882                 double x6 = coord[dim_space*N6];
1883                 double x7 = coord[dim_space*N7];
1884                 double x8 = coord[dim_space*N8];
1885
1886                 double y1 = coord[dim_space*N1+1];
1887                 double y2 = coord[dim_space*N2+1];
1888                 double y3 = coord[dim_space*N3+1];
1889                 double y4 = coord[dim_space*N4+1];
1890                 double y5 = coord[dim_space*N5+1];
1891                 double y6 = coord[dim_space*N6+1];
1892                 double y7 = coord[dim_space*N7+1];
1893                 double y8 = coord[dim_space*N8+1];
1894
1895                 double z1 = coord[dim_space*N1+2];
1896                 double z2 = coord[dim_space*N2+2];
1897                 double z3 = coord[dim_space*N3+2];
1898                 double z4 = coord[dim_space*N4+2];
1899                 double z5 = coord[dim_space*N5+2];
1900                 double z6 = coord[dim_space*N6+2];
1901                 double z7 = coord[dim_space*N7+2];
1902                 double z8 = coord[dim_space*N8+2];
1903
1904                 xbarycenter1 = (x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8)/8.0;
1905                 xbarycenter2 = (y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8)/8.0;
1906                 xbarycenter3 = (z1 + z2 + z3 + z4 + z5 + z6 + z7 + z8)/8.0;
1907
1908                 barycenter[3*index] = xbarycenter1 ;
1909                 barycenter[3*index+1] = xbarycenter2 ;
1910                 barycenter[3*index+2] = xbarycenter3 ;
1911
1912                 index++;
1913               }
1914             break;
1915           }
1916         case MED_TRIA3 : case MED_TRIA6 :
1917           {
1918             for (int tria=0;tria<nb_entity_type;tria++)
1919               {
1920                 int tria_index = (type%100)*tria;
1921
1922                 int N1 = global_connectivity[tria_index]-1;
1923                 int N2 = global_connectivity[tria_index+1]-1;
1924                 int N3 = global_connectivity[tria_index+2]-1;
1925
1926                 double x1 = coord[dim_space*N1];
1927                 double x2 = coord[dim_space*N2];
1928                 double x3 = coord[dim_space*N3];
1929
1930                 double y1 = coord[dim_space*N1+1];
1931                 double y2 = coord[dim_space*N2+1];
1932                 double y3 = coord[dim_space*N3+1];
1933
1934                 xbarycenter1 = (x1 + x2 + x3)/3.0;
1935                 xbarycenter2 = (y1 + y2 + y3)/3.0;
1936
1937                 if (dim_space==2)
1938                   {
1939                     barycenter[2*index] = xbarycenter1 ;
1940                     barycenter[2*index+1] = xbarycenter2 ;
1941                   }
1942                 else
1943                   {
1944                     double z1 =
1945                       coord[dim_space*N1+2];
1946                     double z2 =
1947                       coord[dim_space*N2+2];
1948                     double z3 =
1949                       coord[dim_space*N3+2];
1950
1951                     xbarycenter3 = (z1 + z2 + z3)/3.0;
1952
1953                     barycenter[3*index] = xbarycenter1 ;
1954                     barycenter[3*index+1] = xbarycenter2 ;
1955                     barycenter[3*index+2] = xbarycenter3 ;
1956                   }
1957
1958                 index++;
1959               }
1960             break;
1961           }
1962         case MED_QUAD4 : case MED_QUAD8 :
1963           {
1964             for (int quad=0;quad<nb_entity_type;quad++)
1965               {
1966                 int quad_index = (type%100)*quad;
1967
1968                 int N1 = global_connectivity[quad_index]-1;
1969                 int N2 = global_connectivity[quad_index+1]-1;
1970                 int N3 = global_connectivity[quad_index+2]-1;
1971                 int N4 = global_connectivity[quad_index+3]-1;
1972
1973                 double x1 = coord[dim_space*N1];
1974                 double x2 = coord[dim_space*N2];
1975                 double x3 = coord[dim_space*N3];
1976                 double x4 = coord[dim_space*N4];
1977
1978                 double y1 = coord[dim_space*N1+1];
1979                 double y2 = coord[dim_space*N2+1];
1980                 double y3 = coord[dim_space*N3+1];
1981                 double y4 = coord[dim_space*N4+1];
1982
1983                 xbarycenter1 = (x1 + x2 + x3 + x4)/4.0;
1984                 xbarycenter2 = (y1 + y2 + y3 + y4)/4.0;
1985
1986                 if (dim_space==2)
1987                   {
1988                     barycenter[2*index] = xbarycenter1 ;
1989                     barycenter[2*index+1] = xbarycenter2 ;
1990                   }
1991                 else
1992                   {
1993                     double z1 = coord[dim_space*N1+2];
1994                     double z2 = coord[dim_space*N2+2];
1995                     double z3 = coord[dim_space*N3+2];
1996                     double z4 = coord[dim_space*N4+2];
1997
1998                     xbarycenter3 = (z1 + z2 + z3 + z4)/4.0;
1999
2000                     barycenter[3*index] = xbarycenter1 ;
2001                     barycenter[3*index+1] = xbarycenter2 ;
2002                     barycenter[3*index+2] = xbarycenter3 ;
2003                   }
2004                 index++;
2005               }
2006             break;
2007           }
2008         case MED_SEG2 : case MED_SEG3 :
2009           {
2010             for (int edge=0;edge<nb_entity_type;edge++)
2011               {
2012                 int edge_index = (type%100)*edge;
2013
2014                 int N1 = global_connectivity[edge_index]-1;
2015                 int N2 = global_connectivity[edge_index+1]-1;
2016
2017                 double x1 = coord[dim_space*N1];
2018                 double x2 = coord[dim_space*N2];
2019
2020                 double y1 = coord[dim_space*N1+1];
2021                 double y2 = coord[dim_space*N2+1];
2022
2023                 xbarycenter1 = (x1 + x2)/2.0;
2024                 xbarycenter2 = (y1 + y2)/2.0;
2025
2026                 if (dim_space==2)
2027                   {
2028                     barycenter[2*index] = xbarycenter1 ;
2029                     barycenter[2*index+1] = xbarycenter2 ;
2030                   }
2031                 else
2032                   {
2033                     double z1 = coord[dim_space*N1+2];
2034                     double z2 = coord[dim_space*N2+2];
2035
2036                     xbarycenter3 = (z1 + z2)/2.0;
2037
2038                     barycenter[3*index] = xbarycenter1 ;
2039                     barycenter[3*index+1] = xbarycenter2 ;
2040                     barycenter[3*index+2] = xbarycenter3 ;
2041                   }
2042                 index++;
2043               }
2044             break;
2045           }
2046         default :
2047           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get a barycenter on it (in fact unknown type) !"));
2048           break;
2049         }
2050     }
2051
2052   Barycenter->setValue(MED_FULL_INTERLACE,barycenter);
2053
2054   delete[] barycenter ;
2055
2056   END_OF(LOC);
2057
2058   return Barycenter;
2059 }
2060
2061 //=======================================================================
2062 //function : checkGridFillCoords
2063 //purpose  : if this->_isAGrid, assure that _coordinate is filled
2064 //=======================================================================
2065
2066 inline void MESH::checkGridFillCoords() const
2067 {
2068   if (_isAGrid)
2069     ((GRID *) this)->fillCoordinates();
2070 }
2071
2072 //=======================================================================
2073 //function : checkGridFillConnectivity
2074 //purpose  : if this->_isAGrid, assure that _connectivity is filled
2075 //=======================================================================
2076
2077 inline void MESH::checkGridFillConnectivity() const
2078 {
2079   if (_isAGrid)
2080     ((GRID *) this)->fillConnectivity();
2081 }
2082
2083
2084 void MESH::read(int index)
2085 {
2086   const char * LOC = "MESH::read(int index=0) : ";
2087   BEGIN_OF(LOC);
2088
2089   if (_drivers[index]) {
2090     _drivers[index]->open();
2091     _drivers[index]->read();
2092     _drivers[index]->close();
2093   }
2094   else
2095     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
2096                                      << "The index given is invalid, index must be between  0 and |"
2097                                      << _drivers.size()
2098                                      )
2099                           );
2100   if (_isAGrid)
2101     ((GRID *) this)->fillMeshAfterRead();
2102
2103   END_OF(LOC);
2104 }
2105 //=======================================================================
2106 //function : getSkin
2107 //purpose  :
2108 //=======================================================================
2109
2110 SUPPORT * MESH::getSkin(const SUPPORT * Support3D) throw (MEDEXCEPTION)
2111 {
2112   const char * LOC = "MESH::getSkin : " ;
2113   BEGIN_OF(LOC) ;
2114   // some test :
2115   if (this != Support3D->getMesh())
2116     throw MEDEXCEPTION(STRING(LOC) <<  "no compatibility between *this and SUPPORT::_mesh !");
2117   if (_meshDimension != 3 || Support3D->getEntity() != MED_CELL)
2118       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Defined on 3D cells only"));
2119
2120   // well, all rigth !
2121   SUPPORT * mySupport = new SUPPORT(this,"Skin",MED_FACE);
2122   mySupport->setAll(false);
2123
2124   list<int> myElementsList ;
2125   int i,j, size = 0 ;
2126
2127   calculateConnectivity(MED_FULL_INTERLACE, MED_DESCENDING, MED_CELL);
2128   if (Support3D->isOnAllElements())
2129   {
2130     int * myConnectivityValue = const_cast <int*> (getReverseConnectivity(MED_DESCENDING)) ;
2131     int nbFaces = getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS);
2132     for (i=0, j=1 ; j<=nbFaces; ++j, i += 2)
2133     {
2134       int cellNb1 = myConnectivityValue [i];
2135       int cellNb2 = myConnectivityValue [i+1];
2136       //MESSAGE( " FACE # " << j << " -- Cells: " << cellNb1 << ", " << cellNb2 );
2137       if ((cellNb1 == 0 || cellNb2 == 0) && (cellNb1 + cellNb2 > 0))
2138       {
2139         myElementsList.push_back( j ) ;
2140         size++ ;
2141       }
2142     }
2143   }
2144   else
2145   {
2146     map<int,int> FaceNbEncounterNb;
2147
2148     int * myConnectivityValue = const_cast <int *>
2149       (getConnectivity(MED_FULL_INTERLACE, MED_DESCENDING,
2150                        MED_CELL, MED_ALL_ELEMENTS));
2151     int * myConnectivityIndex = const_cast <int *> (getConnectivityIndex(MED_DESCENDING, MED_CELL));
2152     int * myCellNbs = const_cast <int *> (Support3D->getnumber()->getValue());
2153     int nbCells = Support3D->getnumber()->getLength();
2154     for (i=0; i<nbCells; ++i)
2155     {
2156       int cellNb = myCellNbs [ i ];
2157       int faceFirst = myConnectivityIndex[ cellNb-1 ];
2158       int faceLast  = myConnectivityIndex[ cellNb ];
2159       for (j = faceFirst; j < faceLast; ++j)
2160       {
2161         int faceNb = abs( myConnectivityValue [ j-1 ] );
2162         //MESSAGE( "Cell # " << i << " -- Face: " << faceNb);
2163         if (FaceNbEncounterNb.find( faceNb ) == FaceNbEncounterNb.end())
2164           FaceNbEncounterNb[ faceNb ] = 1;
2165         else
2166           FaceNbEncounterNb[ faceNb ] ++;
2167       }
2168     }
2169     map<int,int>::iterator FaceNbEncounterNbItr;
2170     for (FaceNbEncounterNbItr = FaceNbEncounterNb.begin();
2171          FaceNbEncounterNbItr != FaceNbEncounterNb.end();
2172          FaceNbEncounterNbItr ++)
2173       if ((*FaceNbEncounterNbItr).second == 1)
2174       {
2175         myElementsList.push_back( (*FaceNbEncounterNbItr).first) ;
2176         size++ ;
2177       }
2178   }
2179   // Well, we must know how many geometric type we have found
2180   int * myListArray = new int[size] ;
2181   int id = 0 ;
2182   list<int>::iterator myElementsListIt ;
2183   for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
2184     myListArray[id]=(*myElementsListIt) ;
2185     id ++ ;
2186   }
2187
2188   int numberOfGeometricType ;
2189   medGeometryElement* geometricType ;
2190   int * numberOfGaussPoint ;
2191   int * geometricTypeNumber ;
2192   int * numberOfEntities ;
2193   //  MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY() ;
2194   int * mySkyLineArrayIndex ;
2195
2196   int numberOfType = getNumberOfTypes(MED_FACE) ;
2197   if (numberOfType == 1) { // wonderfull : it's easy !
2198     numberOfGeometricType = 1 ;
2199     geometricType = new medGeometryElement[1] ;
2200     const medGeometryElement *  allType = getTypes(MED_FACE);
2201     geometricType[0] = allType[0] ;
2202     numberOfGaussPoint = new int[1] ;
2203     numberOfGaussPoint[0] = 1 ;
2204     geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
2205     geometricTypeNumber[0] = 0 ;
2206     numberOfEntities = new int[1] ;
2207     numberOfEntities[0] = size ;
2208     mySkyLineArrayIndex = new int[2] ;
2209     mySkyLineArrayIndex[0]=1 ;
2210     mySkyLineArrayIndex[1]=1+size ;
2211   }
2212   else {// hemmm
2213     map<medGeometryElement,int> theType ;
2214     for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
2215       medGeometryElement myType = getElementType(MED_FACE,*myElementsListIt) ;
2216       if (theType.find(myType) != theType.end() )
2217         theType[myType]+=1 ;
2218       else
2219         theType[myType]=1 ;
2220     }
2221     numberOfGeometricType = theType.size() ;
2222     geometricType = new medGeometryElement[numberOfGeometricType] ;
2223     const medGeometryElement *  allType = getTypes(MED_FACE);
2224     numberOfGaussPoint = new int[numberOfGeometricType] ;
2225     geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
2226     numberOfEntities = new int[numberOfGeometricType] ;
2227     mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
2228     int index = 0 ;
2229     mySkyLineArrayIndex[0]=1 ;
2230     map<medGeometryElement,int>::iterator theTypeIt ;
2231     for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
2232       geometricType[index] = (*theTypeIt).first ;
2233       numberOfGaussPoint[index] = 1 ;
2234       geometricTypeNumber[index] = 0 ;
2235       numberOfEntities[index] = (*theTypeIt).second ;
2236       mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfEntities[index] ;
2237       index++ ;
2238     }
2239   }
2240   //  mySkyLineArray->setMEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
2241   MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
2242
2243   mySupport->setNumberOfGeometricType(numberOfGeometricType) ;
2244   mySupport->setGeometricType(geometricType) ;
2245   mySupport->setNumberOfGaussPoint(numberOfGaussPoint) ;
2246   //  mySupport->setGeometricTypeNumber(geometricTypeNumber) ;
2247   mySupport->setNumberOfElements(numberOfEntities) ;
2248   mySupport->setTotalNumberOfElements(size) ;
2249   mySupport->setNumber(mySkyLineArray) ;
2250
2251   delete[] numberOfEntities;
2252   delete[] geometricTypeNumber;
2253   delete[] numberOfGaussPoint;
2254   delete[] geometricType;
2255   delete[] mySkyLineArrayIndex;
2256   delete[] myListArray;
2257 //   delete mySkyLineArray;
2258
2259   END_OF(LOC) ;
2260   return mySupport ;
2261
2262 }