src/MEDMEM/MEDMEM_Mesh.cxx

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