src/MEDMEM/MEDMEM_Mesh.cxx

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