src/MEDMEM/MEDMEM_Mesh.cxx

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