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