src/MEDMEM/MEDMEM_Mesh.cxx

   1 using namespace std;
   2 /*
   3  File Mesh.cxx
   4  $Header$
   5 */
   6
   7 #include <math.h>
   8
   9 #include <list>
  10 #include <map>
  11
  12 #include "MEDMEM_DriversDef.hxx"
  13 #include "MEDMEM_Field.hxx"
  14 #include "MEDMEM_Mesh.hxx"
  15
  16 #include "MEDMEM_Support.hxx"
  17 #include "MEDMEM_Family.hxx"
  18 #include "MEDMEM_Group.hxx"
  19 #include "MEDMEM_Coordinate.hxx"
  20 #include "MEDMEM_Connectivity.hxx"
  21 #include "MEDMEM_CellModel.hxx"
  22 using namespace MEDMEM;
  23 //#include "MEDMEM_Grid.hxx" this inclision should have never be here !!!
  24
  25 //update Families with content list
  26 //int family_count(int family_number, int count, int * entities_number, int * entities_list) ;
  27
  28 // ------- Drivers Management Part
  29
  30 // MESH::INSTANCE_DE<MED_MESH_RDONLY_DRIVER> MESH::inst_med_rdonly ;
  31 //const MESH::INSTANCE * const MESH::instances[] = { &MESH::inst_med_rdonly , &MESH::inst_med_rdwr } ;
  32
  33 // Add a similar line for your personnal driver (step 3)
  34
  35 MESH::INSTANCE_DE<MED_MESH_RDWR_DRIVER>   MESH::inst_med   ;
  36 MESH::INSTANCE_DE<GIBI_MESH_RDWR_DRIVER>  MESH::inst_gibi ;
  37 MESH::INSTANCE_DE<PORFLOW_MESH_RDWR_DRIVER>  MESH::inst_porflow ;
  38 MESH::INSTANCE_DE<VTK_MESH_DRIVER> MESH::inst_vtk;
  39
  40 // Add your own driver in the driver list       (step 4)
  41 // Note the list must be coherent with the driver type list defined in MEDMEM_DRIVER.hxx.
  42 const MESH::INSTANCE * const MESH::instances[] =   {  &MESH::inst_med, &MESH::inst_gibi, &MESH::inst_porflow, &MESH::inst_vtk } ;
  43
  44 /*! Add a %MESH driver of type %driverTypes (MED_DRIVER, ....) associated with file fileName. The meshname used in the file
  45     is  driverName. addDriver returns an integer handler. */
  46 int MESH::addDriver(driverTypes driverType,
  47                     const string & fileName/*="Default File Name.med"*/,const string & driverName/*="Default Mesh Name"*/) {
  48
  49   const char * LOC = "MESH::addDriver(driverTypes driverType, const string & fileName=\"Default File Name.med\",const string & driverName=\"Default Mesh Name\") : ";
  50
  51   GENDRIVER * driver;
  52   int current;
  53   int itDriver = (int) NO_DRIVER;
  54
  55   BEGIN_OF(LOC);
  56
  57   SCRUTE(driverType);
  58
  59   SCRUTE(instances[driverType]);
  60
  61   switch(driverType)
  62     {
  63     case MED_DRIVER : {
  64       itDriver = (int) driverType ;
  65       break ;
  66     }
  67
  68     case GIBI_DRIVER : {
  69       itDriver = (int) driverType ;
  70       break ;
  71     }
  72
  73     case PORFLOW_DRIVER : {
  74       itDriver = (int) driverType ;
  75       break ;
  76     }
  77
  78     case VTK_DRIVER : {
  79       itDriver = 3 ;
  80       break ;
  81     }
  82
  83     case NO_DRIVER : {
  84       throw MED_EXCEPTION (LOCALIZED(STRING(LOC)<< "NO_DRIVER has been specified to the method which is not allowed"));
  85     }
  86     }
  87
  88   if (itDriver == ((int) NO_DRIVER))
  89     throw MED_EXCEPTION (LOCALIZED(STRING(LOC)<< "othe driver than MED_DRIVER GIBI_DRIVER PORFLOW_DRIVER and VT_DRIVER has been specified to the method which is not allowed"));
  90
  91   driver = instances[itDriver]->run(fileName, this) ;
  92
  93   _drivers.push_back(driver);
  94
  95   current = _drivers.size()-1;
  96
  97   _drivers[current]->setMeshName(driverName);
  98
  99   END_OF(LOC);
 100
 101   return current;
 102 }
 103
 104 /*! Add an existing MESH driver. */
 105 int  MESH::addDriver(GENDRIVER & driver) {
 106   const char * LOC = "MESH::addDriver(GENDRIVER &) : ";
 107   BEGIN_OF(LOC);
 108
 109   // A faire : Vérifier que le driver est de type MESH.
 110   GENDRIVER * newDriver = driver.copy() ;
 111
 112   _drivers.push_back(newDriver);
 113   return _drivers.size()-1;
 114
 115   END_OF(LOC);
 116 }
 117
 118 /*! Remove an existing MESH driver. */
 119 void MESH::rmDriver (int index/*=0*/) {
 120   const char * LOC = "MESH::rmDriver (int index=0): ";
 121   BEGIN_OF(LOC);
 122
 123   if ( _drivers[index] ) {
 124     //_drivers.erase(&_drivers[index]);
 125     // why not ????
 126     MESSAGE ("detruire");
 127   }
 128   else
 129     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
 130                                      << "The index given is invalid, index must be between  0 and  |"
 131                                      << _drivers.size()
 132                                      )
 133                           );
 134
 135   END_OF(LOC);
 136
 137 };
 138
 139 // ------ End of Drivers Management Part
 140
 141
 142 void MESH::init() {
 143
 144   const char * LOC = "MESH::init(): ";
 145
 146   BEGIN_OF(LOC);
 147
 148   string        _name = "NOT DEFINED"; // A POSITIONNER EN FCT DES IOS ?
 149
 150   string        _decription = "NOT DEFINED"; // A POSITIONNER EN FCT DES IOS ?
 151
 152   _coordinate   = (COORDINATE   *) NULL;
 153   _connectivity = (CONNECTIVITY *) NULL;
 154
 155   _spaceDimension        =          MED_INVALID; // 0 ?!?
 156   _meshDimension         =          MED_INVALID;
 157   _numberOfNodes         =          MED_INVALID;
 158
 159   _isAGrid = false;
 160
 161   _arePresentOptionnalNodesNumbers = 0;
 162
 163   END_OF(LOC);
 164 };
 165
 166 /*! Create an empty MESH. */
 167 MESH::MESH():_coordinate(NULL),_connectivity(NULL), _isAGrid(false) {
 168   init();
 169 };
 170
 171 MESH::MESH(MESH &m)
 172 {
 173   _name = m._name;
 174   _description = m._description;
 175   _isAGrid = m._isAGrid;
 176
 177   if (m._coordinate != NULL)
 178     _coordinate = new COORDINATE(* m._coordinate);
 179   else
 180     _coordinate = (COORDINATE *) NULL;
 181
 182   if (m._connectivity != NULL)
 183     _connectivity = new CONNECTIVITY(* m._connectivity);
 184   else
 185     _connectivity = (CONNECTIVITY *) NULL;
 186
 187   _spaceDimension = m._spaceDimension;
 188   _meshDimension = m._meshDimension;
 189   _numberOfNodes = m._numberOfNodes;
 190
 191   _familyNode = m._familyNode;
 192   for (int i=0; i<m._familyNode.size(); i++)
 193     {
 194       _familyNode[i] = new FAMILY(* m._familyNode[i]);
 195       _familyNode[i]->setMesh(this);
 196     }
 197
 198   _familyCell = m._familyCell;
 199   for (int i=0; i<m._familyCell.size(); i++)
 200     {
 201       _familyCell[i] = new FAMILY(* m._familyCell[i]);
 202       _familyCell[i]->setMesh(this);
 203     }
 204
 205   _familyFace = m._familyFace;
 206   for (int i=0; i<m._familyFace.size(); i++)
 207     {
 208       _familyFace[i] = new FAMILY(* m._familyFace[i]);
 209       _familyFace[i]->setMesh(this);
 210     }
 211
 212   _familyEdge = m._familyEdge;
 213   for (int i=0; i<m._familyEdge.size(); i++)
 214     {
 215       _familyEdge[i] = new FAMILY(* m._familyEdge[i]);
 216       _familyEdge[i]->setMesh(this);
 217     }
 218
 219   _groupNode = m._groupNode;
 220   for (int i=0; i<m._groupNode.size(); i++)
 221     {
 222       _groupNode[i] = new GROUP(* m._groupNode[i]);
 223       _groupNode[i]->setMesh(this);
 224     }
 225
 226   _groupCell = m._groupCell;
 227   for (int i=0; i<m._groupCell.size(); i++)
 228     {
 229       _groupCell[i] = new GROUP(* m._groupCell[i]);
 230       _groupCell[i]->setMesh(this);
 231     }
 232
 233   _groupFace = m._groupFace;
 234   for (int i=0; i<m._groupFace.size(); i++)
 235     {
 236       _groupFace[i] = new GROUP(* m._groupFace[i]);
 237       _groupFace[i]->setMesh(this);
 238     }
 239
 240   _groupEdge = m._groupEdge;
 241   for (int i=0; i<m._groupEdge.size(); i++)
 242     {
 243       _groupEdge[i] = new GROUP(* m._groupEdge[i]);
 244       _groupEdge[i]->setMesh(this);
 245     }
 246
 247   //_drivers = m._drivers;  //Recopie des drivers?
 248 }
 249
 250 MESH::~MESH() {
 251
 252   MESSAGE("MESH::~MESH() : Destroying the Mesh");
 253   if (_coordinate != ((COORDINATE *) NULL)) delete _coordinate ;
 254   if (_connectivity != ((CONNECTIVITY *) NULL)) delete _connectivity ;
 255   int size ;
 256   size = _familyNode.size() ;
 257   for (int i=0;i<size;i++)
 258     delete _familyNode[i] ;
 259   size = _familyCell.size() ;
 260   for (int i=0;i<size;i++)
 261     delete _familyCell[i] ;
 262   size = _familyFace.size() ;
 263   for (int i=0;i<size;i++)
 264     delete _familyFace[i] ;
 265   size = _familyEdge.size() ;
 266   for (int i=0;i<size;i++)
 267     delete _familyEdge[i] ;
 268   size = _groupNode.size() ;
 269   for (int i=0;i<size;i++)
 270     delete _groupNode[i] ;
 271   size = _groupCell.size() ;
 272   for (int i=0;i<size;i++)
 273     delete _groupCell[i] ;
 274   size = _groupFace.size() ;
 275   for (int i=0;i<size;i++)
 276     delete _groupFace[i] ;
 277   size = _groupEdge.size() ;
 278   for (int i=0;i<size;i++)
 279     delete _groupEdge[i] ;
 280
 281   MESSAGE("In this object MESH there is(are) " << _drivers.size() << " driver(s)");
 282
 283   for (unsigned int index=0; index < _drivers.size(); index++ )
 284     {
 285       SCRUTE(_drivers[index]);
 286       if ( _drivers[index] != NULL) delete _drivers[index];
 287     }
 288
 289 }
 290
 291 MESH & MESH::operator=(const MESH &m)
 292 {
 293   const char * LOC = "MESH & MESH::operator=(const MESH &m) : ";
 294   BEGIN_OF(LOC);
 295
 296   MESSAGE(LOC <<"Not yet implemented, operating on the object " << m);
 297   //  A FAIRE.........
 298
 299   // ATTENTION CET OPERATEUR DE RECOPIE EST DANGEREUX POUR LES
 300   // POINTEURS : ex : nodal_connectivity ???? EXPRES ????
 301 //    _drivers = m._drivers;
 302
 303 //    space_dimension=m.space_dimension;
 304 //    mesh_dimension=m.mesh_dimension;
 305
 306 //    nodes_count=m.nodes_count;
 307 //    coordinates=m.coordinates;
 308 //    coordinates_name=m.coordinates_name ;
 309 //    coordinates_unit=m.coordinates_unit ;
 310 //    nodes_numbers=m.nodes_numbers ;
 311
 312 //    cells_types_count=m.cells_types_count;
 313 //    cells_type=m.cells_type;
 314 //    cells_count=m.cells_count;
 315 //    nodal_connectivity=m.nodal_connectivity;
 316
 317 //    nodes_families_count=m.nodes_families_count;
 318 //    nodes_Families=m.nodes_Families;
 319
 320 //    cells_families_count=m.cells_families_count;
 321 //    cells_Families=m.cells_Families;
 322
 323 //    maximum_cell_number_by_node = m.maximum_cell_number_by_node;
 324 //    if (maximum_cell_number_by_node > 0)
 325 //      {
 326 //        reverse_nodal_connectivity = m.reverse_nodal_connectivity;
 327 //        reverse_nodal_connectivity_index = m.reverse_nodal_connectivity_index ;
 328 //      }
 329   END_OF(LOC);
 330
 331   return *this;
 332 }
 333
 334 /*! Create a %MESH object using a %MESH driver of type %driverTypes (MED_DRIVER, ....) associated with file fileName.
 335   The meshname driverName must already exists in the file.*/
 336 MESH::MESH(driverTypes driverType, const string &  fileName/*=""*/, const string &  driverName/*=""*/) throw (MEDEXCEPTION)
 337 {
 338   const char * LOC ="MESH::MESH(driverTypes driverType, const string &  fileName="", const string &  driverName="") : ";
 339   int current;
 340
 341   BEGIN_OF(LOC);
 342
 343   init();
 344
 345   switch(driverType)
 346     {
 347     case MED_DRIVER :
 348       {
 349         MED_MESH_RDONLY_DRIVER myDriver(fileName,this);
 350         myDriver.setMeshName(driverName);
 351         current = addDriver(myDriver);
 352         break;
 353       }
 354     case GIBI_DRIVER :
 355       {
 356         GIBI_MESH_RDONLY_DRIVER myDriver(fileName,this);
 357         current = addDriver(myDriver);
 358         break;
 359       }
 360     case PORFLOW_DRIVER :
 361       {
 362         PORFLOW_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 & MEDMEM::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(((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3))*
1587                               ((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3)) +
1588                               ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3))*
1589                               ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3)) +
1590                               ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))*
1591                               ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))))/2.0;
1592
1593 //                       sqrt(((y3-y2)*(z4-z3) - (y4-y3)*(z3-z2))*
1594 //                            ((y3-y2)*(z4-z3) - (y4-y3)*(z3-z2)) +
1595 //                            ((x4-x3)*(z2-z2) - (x3-x2)*(z4-z3))*
1596 //                            ((x4-x3)*(z4-z2) - (x3-x2)*(z4-z3)) +
1597 //                            ((x3-x2)*(y4-y3) - (x4-x3)*(y3-y2))*
1598 //                            ((x3-x2)*(y4-y3) - (x4-x3)*(y3-y2))))/2.0;
1599
1600                 xnormal1 = xnormal1*xarea;
1601                 xnormal2 = xnormal2*xarea;
1602                 xnormal3 = xnormal3*xarea;
1603
1604                 normal[3*index] = xnormal1 ;
1605                 normal[3*index+1] = xnormal2 ;
1606                 normal[3*index+2] = xnormal3 ;
1607
1608                 index++;
1609               }
1610             break;
1611           }
1612         case MED_SEG2 : case MED_SEG3 :
1613           {
1614             for (int edge=0;edge<nb_entity_type;edge++)
1615               {
1616                 int edge_index = (type%100)*edge;
1617
1618                 int N1 = global_connectivity[edge_index]-1;
1619                 int N2 = global_connectivity[edge_index+1]-1;
1620
1621                 double x1 = coord[dim_space*N1];
1622                 double x2 = coord[dim_space*N2];
1623
1624                 double y1 = coord[dim_space*N1+1];
1625                 double y2 = coord[dim_space*N2+1];
1626
1627                 xnormal1 = -(y2-y1);
1628                 xnormal2 = x2-x1;
1629
1630                 normal[2*index] = xnormal1 ;
1631                 normal[2*index+1] = xnormal2 ;
1632
1633                 index++;
1634               }
1635             break;
1636           }
1637         default :
1638           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Normals on it !"));
1639           break;
1640         }
1641
1642       if (!onAll) delete [] global_connectivity ;
1643     }
1644
1645   Normal->setValue(MED_FULL_INTERLACE,normal);
1646   delete[] normal ;
1647
1648   END_OF(LOC);
1649
1650   return Normal;
1651 }
1652
1653 FIELD<double>* MESH::getBarycenter(const SUPPORT * Support) const throw (MEDEXCEPTION)
1654 {
1655   const char * LOC = "MESH::getBarycenter(SUPPORT*) : ";
1656   BEGIN_OF(LOC);
1657
1658   // Make sure that the MESH class is the same as the MESH class attribut
1659   // in the class Support
1660   MESH* myMesh = Support->getMesh();
1661   if (this != myMesh)
1662     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1663
1664   int dim_space = getSpaceDimension();
1665   medEntityMesh support_entity = Support->getEntity();
1666   bool onAll = Support->isOnAllElements();
1667
1668   int nb_type, length_values;
1669   const medGeometryElement* types;
1670   int nb_entity_type;
1671   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1672   const int* global_connectivity;
1673
1674 //    if (onAll)
1675 //      {
1676 //        nb_type = myMesh->getNumberOfTypes(support_entity);
1677 //        length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1678 //        types = getTypes(support_entity);
1679 //      }
1680 //    else
1681     {
1682       nb_type = Support->getNumberOfTypes();
1683       length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1684       types = Support->getTypes();
1685     }
1686
1687   int index;
1688   FIELD<double>* Barycenter;
1689
1690   Barycenter = new FIELD<double>::FIELD(Support,dim_space);
1691   Barycenter->setName("BARYCENTER");
1692   Barycenter->setDescription("cells or faces barycenter");
1693
1694   for (int k=0;k<dim_space;k++) {
1695     int kp1 = k + 1;
1696     Barycenter->setComponentName(kp1,myMesh->getCoordinatesNames()[k]);
1697     Barycenter->setComponentDescription(kp1,"desc-comp");
1698     Barycenter->setMEDComponentUnit(kp1,myMesh->getCoordinatesUnits()[k]);
1699   }
1700
1701   Barycenter->setValueType(MED_REEL64);
1702
1703   Barycenter->setIterationNumber(0);
1704   Barycenter->setOrderNumber(0);
1705   Barycenter->setTime(0.0);
1706
1707   double *barycenter = new double [dim_space*length_values];
1708   //  double *barycenter = Barycenter->getValue(MED_FULL_INTERLACE);
1709
1710   const double * coord = getCoordinates(MED_FULL_INTERLACE);
1711   index = 0;
1712
1713   for (int i=0;i<nb_type;i++)
1714     {
1715       medGeometryElement type = types[i] ;
1716       double xbarycenter1, xbarycenter2, xbarycenter3;
1717
1718       if (onAll)
1719         {
1720           nb_entity_type = getNumberOfElements(support_entity,type);
1721           global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1722         }
1723       else
1724         {
1725           // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all !"));
1726           nb_entity_type = Support->getNumberOfElements(type);
1727
1728           const int * supp_number = Support->getNumber(type);
1729           const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1730           const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1731           int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1732
1733           for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1734             for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1735               global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1736             }
1737           }
1738           global_connectivity = global_connectivity_tmp;
1739         }
1740
1741       switch (type)
1742         {
1743         case MED_TETRA4 : case MED_TETRA10 :
1744           {
1745             for (int tetra =0;tetra<nb_entity_type;tetra++)
1746               {
1747                 int tetra_index = (type%100)*tetra;
1748
1749                 int N1 = global_connectivity[tetra_index]-1;
1750                 int N2 = global_connectivity[tetra_index+1]-1;
1751                 int N3 = global_connectivity[tetra_index+2]-1;
1752                 int N4 = global_connectivity[tetra_index+3]-1;
1753
1754                 double x1 = coord[dim_space*N1];
1755                 double x2 = coord[dim_space*N2];
1756                 double x3 = coord[dim_space*N3];
1757                 double x4 = coord[dim_space*N4];
1758
1759                 double y1 = coord[dim_space*N1+1];
1760                 double y2 = coord[dim_space*N2+1];
1761                 double y3 = coord[dim_space*N3+1];
1762                 double y4 = coord[dim_space*N4+1];
1763
1764                 double z1 = coord[dim_space*N1+2];
1765                 double z2 = coord[dim_space*N2+2];
1766                 double z3 = coord[dim_space*N3+2];
1767                 double z4 = coord[dim_space*N4+2];
1768
1769                 xbarycenter1 = (x1 + x2 + x3 + x4)/4.0;
1770                 xbarycenter2 = (y1 + y2 + y3 + y4)/4.0;
1771                 xbarycenter3 = (z1 + z2 + z3 + z4)/4.0;
1772                 barycenter[3*index] = xbarycenter1 ;
1773                 barycenter[3*index+1] = xbarycenter2 ;
1774                 barycenter[3*index+2] = xbarycenter3 ;
1775                 index++;
1776               }
1777             break;
1778           }
1779         case MED_PYRA5 : case MED_PYRA13 :
1780           {
1781             for (int pyra=0;pyra<nb_entity_type;pyra++)
1782               {
1783                 int pyra_index = (type%100)*pyra;
1784
1785                 int N1 = global_connectivity[pyra_index]-1;
1786                 int N2 = global_connectivity[pyra_index+1]-1;
1787                 int N3 = global_connectivity[pyra_index+2]-1;
1788                 int N4 = global_connectivity[pyra_index+3]-1;
1789                 int N5 = global_connectivity[pyra_index+4]-1;
1790
1791                 double x1 = coord[dim_space*N1];
1792                 double x2 = coord[dim_space*N2];
1793                 double x3 = coord[dim_space*N3];
1794                 double x4 = coord[dim_space*N4];
1795                 double x5 = coord[dim_space*N5];
1796
1797                 double y1 = coord[dim_space*N1+1];
1798                 double y2 = coord[dim_space*N2+1];
1799                 double y3 = coord[dim_space*N3+1];
1800                 double y4 = coord[dim_space*N4+1];
1801                 double y5 = coord[dim_space*N5+1];
1802
1803                 double z1 = coord[dim_space*N1+2];
1804                 double z2 = coord[dim_space*N2+2];
1805                 double z3 = coord[dim_space*N3+2];
1806                 double z4 = coord[dim_space*N4+2];
1807                 double z5 = coord[dim_space*N5+2];
1808
1809                 xbarycenter1 = (x1 + x2 + x3 + x4 + x5)/5.0;
1810                 xbarycenter2 = (y1 + y2 + y3 + y4 + y5)/5.0;
1811                 xbarycenter3 = (z1 + z2 + z3 + z4 + z5)/5.0;
1812                 barycenter[3*index] = xbarycenter1 ;
1813                 barycenter[3*index+1] = xbarycenter2 ;
1814                 barycenter[3*index+2] = xbarycenter3 ;
1815                 index++;
1816               }
1817             break;
1818           }
1819         case MED_PENTA6 : case MED_PENTA15 :
1820           {
1821             for (int penta=0;penta<nb_entity_type;penta++)
1822               {
1823                 int penta_index = (type%100)*penta;
1824
1825                 int N1 = global_connectivity[penta_index]-1;
1826                 int N2 = global_connectivity[penta_index+1]-1;
1827                 int N3 = global_connectivity[penta_index+2]-1;
1828                 int N4 = global_connectivity[penta_index+3]-1;
1829                 int N5 = global_connectivity[penta_index+4]-1;
1830                 int N6 = global_connectivity[penta_index+5]-1;
1831
1832                 double x1 = coord[dim_space*N1];
1833                 double x2 = coord[dim_space*N2];
1834                 double x3 = coord[dim_space*N3];
1835                 double x4 = coord[dim_space*N4];
1836                 double x5 = coord[dim_space*N5];
1837                 double x6 = coord[dim_space*N6];
1838
1839                 double y1 = coord[dim_space*N1+1];
1840                 double y2 = coord[dim_space*N2+1];
1841                 double y3 = coord[dim_space*N3+1];
1842                 double y4 = coord[dim_space*N4+1];
1843                 double y5 = coord[dim_space*N5+1];
1844                 double y6 = coord[dim_space*N6+1];
1845
1846                 double z1 = coord[dim_space*N1+2];
1847                 double z2 = coord[dim_space*N2+2];
1848                 double z3 = coord[dim_space*N3+2];
1849                 double z4 = coord[dim_space*N4+2];
1850                 double z5 = coord[dim_space*N5+2];
1851                 double z6 = coord[dim_space*N6+2];
1852
1853                 xbarycenter1 = (x1 + x2 + x3 + x4 + x5 + x6)/6.0;
1854                 xbarycenter2 = (y1 + y2 + y3 + y4 + y5 + y6)/6.0;
1855                 xbarycenter3 = (z1 + z2 + z3 + z4 + z5 + z6)/6.0;
1856                 barycenter[3*index] = xbarycenter1 ;
1857                 barycenter[3*index+1] = xbarycenter2 ;
1858                 barycenter[3*index+2] = xbarycenter3 ;
1859                 index++;
1860               }
1861             break;
1862           }
1863         case MED_HEXA8 : case MED_HEXA20 :
1864           {
1865             for (int hexa=0;hexa<nb_entity_type;hexa++)
1866               {
1867                 int hexa_index = (type%100)*hexa;
1868
1869                 int N1 = global_connectivity[hexa_index]-1;
1870                 int N2 = global_connectivity[hexa_index+1]-1;
1871                 int N3 = global_connectivity[hexa_index+2]-1;
1872                 int N4 = global_connectivity[hexa_index+3]-1;
1873                 int N5 = global_connectivity[hexa_index+4]-1;
1874                 int N6 = global_connectivity[hexa_index+5]-1;
1875                 int N7 = global_connectivity[hexa_index+6]-1;
1876                 int N8 = global_connectivity[hexa_index+7]-1;
1877
1878                 double x1 = coord[dim_space*N1];
1879                 double x2 = coord[dim_space*N2];
1880                 double x3 = coord[dim_space*N3];
1881                 double x4 = coord[dim_space*N4];
1882                 double x5 = coord[dim_space*N5];
1883                 double x6 = coord[dim_space*N6];
1884                 double x7 = coord[dim_space*N7];
1885                 double x8 = coord[dim_space*N8];
1886
1887                 double y1 = coord[dim_space*N1+1];
1888                 double y2 = coord[dim_space*N2+1];
1889                 double y3 = coord[dim_space*N3+1];
1890                 double y4 = coord[dim_space*N4+1];
1891                 double y5 = coord[dim_space*N5+1];
1892                 double y6 = coord[dim_space*N6+1];
1893                 double y7 = coord[dim_space*N7+1];
1894                 double y8 = coord[dim_space*N8+1];
1895
1896                 double z1 = coord[dim_space*N1+2];
1897                 double z2 = coord[dim_space*N2+2];
1898                 double z3 = coord[dim_space*N3+2];
1899                 double z4 = coord[dim_space*N4+2];
1900                 double z5 = coord[dim_space*N5+2];
1901                 double z6 = coord[dim_space*N6+2];
1902                 double z7 = coord[dim_space*N7+2];
1903                 double z8 = coord[dim_space*N8+2];
1904
1905                 xbarycenter1 = (x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8)/8.0;
1906                 xbarycenter2 = (y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8)/8.0;
1907                 xbarycenter3 = (z1 + z2 + z3 + z4 + z5 + z6 + z7 + z8)/8.0;
1908
1909                 barycenter[3*index] = xbarycenter1 ;
1910                 barycenter[3*index+1] = xbarycenter2 ;
1911                 barycenter[3*index+2] = xbarycenter3 ;
1912
1913                 index++;
1914               }
1915             break;
1916           }
1917         case MED_TRIA3 : case MED_TRIA6 :
1918           {
1919             for (int tria=0;tria<nb_entity_type;tria++)
1920               {
1921                 int tria_index = (type%100)*tria;
1922
1923                 int N1 = global_connectivity[tria_index]-1;
1924                 int N2 = global_connectivity[tria_index+1]-1;
1925                 int N3 = global_connectivity[tria_index+2]-1;
1926
1927                 double x1 = coord[dim_space*N1];
1928                 double x2 = coord[dim_space*N2];
1929                 double x3 = coord[dim_space*N3];
1930
1931                 double y1 = coord[dim_space*N1+1];
1932                 double y2 = coord[dim_space*N2+1];
1933                 double y3 = coord[dim_space*N3+1];
1934
1935                 xbarycenter1 = (x1 + x2 + x3)/3.0;
1936                 xbarycenter2 = (y1 + y2 + y3)/3.0;
1937
1938                 if (dim_space==2)
1939                   {
1940                     barycenter[2*index] = xbarycenter1 ;
1941                     barycenter[2*index+1] = xbarycenter2 ;
1942                   }
1943                 else
1944                   {
1945                     double z1 =
1946                       coord[dim_space*N1+2];
1947                     double z2 =
1948                       coord[dim_space*N2+2];
1949                     double z3 =
1950                       coord[dim_space*N3+2];
1951
1952                     xbarycenter3 = (z1 + z2 + z3)/3.0;
1953
1954                     barycenter[3*index] = xbarycenter1 ;
1955                     barycenter[3*index+1] = xbarycenter2 ;
1956                     barycenter[3*index+2] = xbarycenter3 ;
1957                   }
1958
1959                 index++;
1960               }
1961             break;
1962           }
1963         case MED_QUAD4 : case MED_QUAD8 :
1964           {
1965             for (int quad=0;quad<nb_entity_type;quad++)
1966               {
1967                 int quad_index = (type%100)*quad;
1968
1969                 int N1 = global_connectivity[quad_index]-1;
1970                 int N2 = global_connectivity[quad_index+1]-1;
1971                 int N3 = global_connectivity[quad_index+2]-1;
1972                 int N4 = global_connectivity[quad_index+3]-1;
1973
1974                 double x1 = coord[dim_space*N1];
1975                 double x2 = coord[dim_space*N2];
1976                 double x3 = coord[dim_space*N3];
1977                 double x4 = coord[dim_space*N4];
1978
1979                 double y1 = coord[dim_space*N1+1];
1980                 double y2 = coord[dim_space*N2+1];
1981                 double y3 = coord[dim_space*N3+1];
1982                 double y4 = coord[dim_space*N4+1];
1983
1984                 xbarycenter1 = (x1 + x2 + x3 + x4)/4.0;
1985                 xbarycenter2 = (y1 + y2 + y3 + y4)/4.0;
1986
1987                 if (dim_space==2)
1988                   {
1989                     barycenter[2*index] = xbarycenter1 ;
1990                     barycenter[2*index+1] = xbarycenter2 ;
1991                   }
1992                 else
1993                   {
1994                     double z1 = coord[dim_space*N1+2];
1995                     double z2 = coord[dim_space*N2+2];
1996                     double z3 = coord[dim_space*N3+2];
1997                     double z4 = coord[dim_space*N4+2];
1998
1999                     xbarycenter3 = (z1 + z2 + z3 + z4)/4.0;
2000
2001                     barycenter[3*index] = xbarycenter1 ;
2002                     barycenter[3*index+1] = xbarycenter2 ;
2003                     barycenter[3*index+2] = xbarycenter3 ;
2004                   }
2005                 index++;
2006               }
2007             break;
2008           }
2009         case MED_SEG2 : case MED_SEG3 :
2010           {
2011             for (int edge=0;edge<nb_entity_type;edge++)
2012               {
2013                 int edge_index = (type%100)*edge;
2014
2015                 int N1 = global_connectivity[edge_index]-1;
2016                 int N2 = global_connectivity[edge_index+1]-1;
2017
2018                 double x1 = coord[dim_space*N1];
2019                 double x2 = coord[dim_space*N2];
2020
2021                 double y1 = coord[dim_space*N1+1];
2022                 double y2 = coord[dim_space*N2+1];
2023
2024                 xbarycenter1 = (x1 + x2)/2.0;
2025                 xbarycenter2 = (y1 + y2)/2.0;
2026
2027                 if (dim_space==2)
2028                   {
2029                     barycenter[2*index] = xbarycenter1 ;
2030                     barycenter[2*index+1] = xbarycenter2 ;
2031                   }
2032                 else
2033                   {
2034                     double z1 = coord[dim_space*N1+2];
2035                     double z2 = coord[dim_space*N2+2];
2036
2037                     xbarycenter3 = (z1 + z2)/2.0;
2038
2039                     barycenter[3*index] = xbarycenter1 ;
2040                     barycenter[3*index+1] = xbarycenter2 ;
2041                     barycenter[3*index+2] = xbarycenter3 ;
2042                   }
2043                 index++;
2044               }
2045             break;
2046           }
2047         default :
2048           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get a barycenter on it (in fact unknown type) !"));
2049           break;
2050         }
2051
2052       if (!onAll) delete [] global_connectivity ;
2053     }
2054
2055   Barycenter->setValue(MED_FULL_INTERLACE,barycenter);
2056
2057   delete[] barycenter ;
2058
2059   END_OF(LOC);
2060
2061   return Barycenter;
2062 }
2063
2064 //=======================================================================
2065 //function : checkGridFillCoords
2066 //purpose  : if this->_isAGrid, assure that _coordinate is filled
2067 //=======================================================================
2068
2069 // inline void MESH::checkGridFillCoords() const
2070 // {
2071 //   if (_isAGrid)
2072 //     ((GRID *) this)->fillCoordinates();
2073 // }
2074
2075 //=======================================================================
2076 //function : checkGridFillConnectivity
2077 //purpose  : if this->_isAGrid, assure that _connectivity is filled
2078 //=======================================================================
2079
2080 // inline void MESH::checkGridFillConnectivity() const
2081 // {
2082 //   if (_isAGrid)
2083 //     ((GRID *) this)->fillConnectivity();
2084 // }
2085
2086 bool MESH::isEmpty() const
2087 {
2088   bool notempty = _name != "" || _description != "" ||
2089                   _coordinate != NULL || _connectivity != NULL ||
2090                   _spaceDimension !=MED_INVALID ||
2091                   _meshDimension !=MED_INVALID ||
2092                   _numberOfNodes !=MED_INVALID || _groupNode.size() != 0 ||
2093                   _familyNode.size() != 0 || _groupCell.size() != 0 ||
2094                   _familyCell.size() != 0 || _groupFace.size() != 0 ||
2095                   _familyFace.size() != 0 || _groupEdge.size() != 0 ||
2096                   _familyEdge.size() != 0 || _isAGrid != 0 ;
2097   return !notempty;
2098 }
2099
2100 void MESH::read(int index)
2101 {
2102   const char * LOC = "MESH::read(int index=0) : ";
2103   BEGIN_OF(LOC);
2104
2105   if (_drivers[index]) {
2106     _drivers[index]->open();
2107     _drivers[index]->read();
2108     _drivers[index]->close();
2109   }
2110   else
2111     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
2112                                      << "The index given is invalid, index must be between  0 and |"
2113                                      << _drivers.size()
2114                                      )
2115                           );
2116 //   if (_isAGrid)
2117 //     ((GRID *) this)->fillMeshAfterRead();
2118
2119   END_OF(LOC);
2120 }
2121 //=======================================================================
2122 //function : getSkin
2123 //purpose  :
2124 //=======================================================================
2125
2126 SUPPORT * MESH::getSkin(const SUPPORT * Support3D) throw (MEDEXCEPTION)
2127 {
2128   const char * LOC = "MESH::getSkin : " ;
2129   BEGIN_OF(LOC) ;
2130   // some test :
2131   if (this != Support3D->getMesh())
2132     throw MEDEXCEPTION(STRING(LOC) <<  "no compatibility between *this and SUPPORT::_mesh !");
2133   if (_meshDimension != 3 || Support3D->getEntity() != MED_CELL)
2134       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Defined on 3D cells only"));
2135
2136   // well, all rigth !
2137   SUPPORT * mySupport = new SUPPORT(this,"Skin",MED_FACE);
2138   mySupport->setAll(false);
2139
2140   list<int> myElementsList ;
2141   int i,j, size = 0 ;
2142
2143   calculateConnectivity(MED_FULL_INTERLACE, MED_DESCENDING, MED_CELL);
2144   if (Support3D->isOnAllElements())
2145   {
2146     int * myConnectivityValue = const_cast <int*> (getReverseConnectivity(MED_DESCENDING)) ;
2147     int nbFaces = getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS);
2148     for (i=0, j=1 ; j<=nbFaces; ++j, i += 2)
2149     {
2150       int cellNb1 = myConnectivityValue [i];
2151       int cellNb2 = myConnectivityValue [i+1];
2152       //MESSAGE( " FACE # " << j << " -- Cells: " << cellNb1 << ", " << cellNb2 );
2153       if ((cellNb1 == 0 || cellNb2 == 0) && (cellNb1 + cellNb2 > 0))
2154       {
2155         myElementsList.push_back( j ) ;
2156         size++ ;
2157       }
2158     }
2159   }
2160   else
2161   {
2162     map<int,int> FaceNbEncounterNb;
2163
2164     int * myConnectivityValue = const_cast <int *>
2165       (getConnectivity(MED_FULL_INTERLACE, MED_DESCENDING,
2166                        MED_CELL, MED_ALL_ELEMENTS));
2167     int * myConnectivityIndex = const_cast <int *> (getConnectivityIndex(MED_DESCENDING, MED_CELL));
2168     int * myCellNbs = const_cast <int *> (Support3D->getnumber()->getValue());
2169     int nbCells = Support3D->getnumber()->getLength();
2170     for (i=0; i<nbCells; ++i)
2171     {
2172       int cellNb = myCellNbs [ i ];
2173       int faceFirst = myConnectivityIndex[ cellNb-1 ];
2174       int faceLast  = myConnectivityIndex[ cellNb ];
2175       for (j = faceFirst; j < faceLast; ++j)
2176       {
2177         int faceNb = abs( myConnectivityValue [ j-1 ] );
2178         //MESSAGE( "Cell # " << i << " -- Face: " << faceNb);
2179         if (FaceNbEncounterNb.find( faceNb ) == FaceNbEncounterNb.end())
2180           FaceNbEncounterNb[ faceNb ] = 1;
2181         else
2182           FaceNbEncounterNb[ faceNb ] ++;
2183       }
2184     }
2185     map<int,int>::iterator FaceNbEncounterNbItr;
2186     for (FaceNbEncounterNbItr = FaceNbEncounterNb.begin();
2187          FaceNbEncounterNbItr != FaceNbEncounterNb.end();
2188          FaceNbEncounterNbItr ++)
2189       if ((*FaceNbEncounterNbItr).second == 1)
2190       {
2191         myElementsList.push_back( (*FaceNbEncounterNbItr).first) ;
2192         size++ ;
2193       }
2194   }
2195   // Well, we must know how many geometric type we have found
2196   int * myListArray = new int[size] ;
2197   int id = 0 ;
2198   list<int>::iterator myElementsListIt ;
2199   for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
2200     myListArray[id]=(*myElementsListIt) ;
2201     id ++ ;
2202   }
2203
2204   int numberOfGeometricType ;
2205   medGeometryElement* geometricType ;
2206   int * numberOfGaussPoint ;
2207   int * geometricTypeNumber ;
2208   int * numberOfEntities ;
2209   //  MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY() ;
2210   int * mySkyLineArrayIndex ;
2211
2212   int numberOfType = getNumberOfTypes(MED_FACE) ;
2213   if (numberOfType == 1) { // wonderfull : it's easy !
2214     numberOfGeometricType = 1 ;
2215     geometricType = new medGeometryElement[1] ;
2216     const medGeometryElement *  allType = getTypes(MED_FACE);
2217     geometricType[0] = allType[0] ;
2218     numberOfGaussPoint = new int[1] ;
2219     numberOfGaussPoint[0] = 1 ;
2220     geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
2221     geometricTypeNumber[0] = 0 ;
2222     numberOfEntities = new int[1] ;
2223     numberOfEntities[0] = size ;
2224     mySkyLineArrayIndex = new int[2] ;
2225     mySkyLineArrayIndex[0]=1 ;
2226     mySkyLineArrayIndex[1]=1+size ;
2227   }
2228   else {// hemmm
2229     map<medGeometryElement,int> theType ;
2230     for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
2231       medGeometryElement myType = getElementType(MED_FACE,*myElementsListIt) ;
2232       if (theType.find(myType) != theType.end() )
2233         theType[myType]+=1 ;
2234       else
2235         theType[myType]=1 ;
2236     }
2237     numberOfGeometricType = theType.size() ;
2238     geometricType = new medGeometryElement[numberOfGeometricType] ;
2239     //const medGeometryElement *  allType = getTypes(MED_FACE); !! UNUSED VARIABLE !!
2240     numberOfGaussPoint = new int[numberOfGeometricType] ;
2241     geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
2242     numberOfEntities = new int[numberOfGeometricType] ;
2243     mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
2244     int index = 0 ;
2245     mySkyLineArrayIndex[0]=1 ;
2246     map<medGeometryElement,int>::iterator theTypeIt ;
2247     for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
2248       geometricType[index] = (*theTypeIt).first ;
2249       numberOfGaussPoint[index] = 1 ;
2250       geometricTypeNumber[index] = 0 ;
2251       numberOfEntities[index] = (*theTypeIt).second ;
2252       mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfEntities[index] ;
2253       index++ ;
2254     }
2255   }
2256   //  mySkyLineArray->setMEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
2257   MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
2258
2259   mySupport->setNumberOfGeometricType(numberOfGeometricType) ;
2260   mySupport->setGeometricType(geometricType) ;
2261   mySupport->setNumberOfGaussPoint(numberOfGaussPoint) ;
2262   //  mySupport->setGeometricTypeNumber(geometricTypeNumber) ;
2263   mySupport->setNumberOfElements(numberOfEntities) ;
2264   mySupport->setTotalNumberOfElements(size) ;
2265   mySupport->setNumber(mySkyLineArray) ;
2266
2267   delete[] numberOfEntities;
2268   delete[] geometricTypeNumber;
2269   delete[] numberOfGaussPoint;
2270   delete[] geometricType;
2271   delete[] mySkyLineArrayIndex;
2272   delete[] myListArray;
2273 //   delete mySkyLineArray;
2274
2275   END_OF(LOC) ;
2276   return mySupport ;
2277
2278 }
2279
2280 /*!
2281   return a SUPPORT pointer on the union of all SUPPORTs in Supports.
2282   You should delete this pointer after use to avoid memory leaks.
2283 */
2284 SUPPORT * MESH::mergeSupports(const vector<SUPPORT *> Supports) const throw (MEDEXCEPTION)
2285 {
2286   const char * LOC = "MESH:::mergeSupports(const vector<SUPPORT *> ) : " ;
2287   BEGIN_OF(LOC) ;
2288
2289   SUPPORT * returnedSupport;
2290   string returnedSupportName;
2291   string returnedSupportDescription;
2292   char * returnedSupportNameChar;
2293   char * returnedSupportDescriptionChar;
2294   int size = Supports.size();
2295
2296   if (size == 1)
2297     {
2298       MESSAGE(LOC <<" there is only one SUPPORT in the argument list, the method return a copy of this object !");
2299       SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2300
2301       returnedSupport = new SUPPORT(*obj);
2302
2303       int lenName = strlen((Supports[0]->getName()).c_str()) + 8 + 1;
2304       int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 8 + 1;
2305
2306       returnedSupportNameChar = new char[lenName];
2307       returnedSupportDescriptionChar = new char[lenDescription];
2308
2309       returnedSupportNameChar = strcpy(returnedSupportNameChar,"Copy of ");
2310       returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
2311       returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Copy of ");
2312       returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2313                                               (Supports[0]->getDescription()).c_str());
2314
2315       returnedSupportName = string(returnedSupportNameChar);
2316       returnedSupportDescription = string(returnedSupportDescriptionChar);
2317
2318       returnedSupport->setName(returnedSupportName);
2319       returnedSupport->setDescription(returnedSupportDescription);
2320
2321       delete [] returnedSupportNameChar;
2322       delete [] returnedSupportDescriptionChar;
2323     }
2324   else
2325     {
2326       SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2327       returnedSupport = new SUPPORT(*obj);
2328
2329       int lenName = strlen((Supports[0]->getName()).c_str()) + 9 + 1;
2330       int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 9 + 1;
2331
2332       for(int i = 1;i<size;i++)
2333         {
2334           obj = const_cast <SUPPORT *> (Supports[i]);
2335           returnedSupport->blending(obj);
2336
2337           if (i == (size-1))
2338             {
2339               lenName = lenName + 5 + strlen((Supports[i]->getName()).c_str());
2340               lenDescription = lenDescription + 5 +
2341                 strlen((Supports[i]->getDescription()).c_str());
2342             }
2343           else
2344             {
2345               lenName = lenName + 2 + strlen((Supports[i]->getName()).c_str());
2346               lenDescription = lenDescription + 2 +
2347                 strlen((Supports[i]->getDescription()).c_str());
2348             }
2349         }
2350
2351       returnedSupportNameChar = new char[lenName];
2352       returnedSupportDescriptionChar = new char[lenDescription];
2353
2354       returnedSupportNameChar = strcpy(returnedSupportNameChar,"Merge of ");
2355       returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Merge of ");
2356
2357       returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
2358       returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2359                                               (Supports[0]->getDescription()).c_str());
2360
2361       for(int i = 1;i<size;i++)
2362         {
2363           if (i == (size-1))
2364             {
2365               returnedSupportNameChar = strcat(returnedSupportNameChar," and ");
2366               returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar," and ");
2367
2368               returnedSupportNameChar = strcat(returnedSupportNameChar,
2369                                                (Supports[i]->getName()).c_str());
2370               returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2371                                                       (Supports[i]->getDescription()).c_str());
2372             }
2373           else
2374             {
2375               returnedSupportNameChar = strcat(returnedSupportNameChar,", ");
2376               returnedSupportNameChar = strcat(returnedSupportNameChar,
2377                                                (Supports[i]->getName()).c_str());
2378
2379               returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,", ");
2380               returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2381                                                       (Supports[i]->getDescription()).c_str());
2382             }
2383         }
2384
2385       returnedSupportName = string(returnedSupportNameChar);
2386       returnedSupport->setName(returnedSupportName);
2387
2388       returnedSupportDescription = string(returnedSupportDescriptionChar);
2389       returnedSupport->setDescription(returnedSupportDescription);
2390
2391       delete [] returnedSupportNameChar;
2392       delete [] returnedSupportDescriptionChar;
2393     }
2394
2395   END_OF(LOC) ;
2396   return returnedSupport;
2397 }
2398
2399 /*!
2400   return a SUPPORT pointer on the intersection of all SUPPORTs in Supports.
2401   The (SUPPORT *) NULL pointer is returned if the intersection is empty.
2402   You should delete this pointer after use to avois memory leaks.
2403 */
2404 SUPPORT * MESH::intersectSupports(const vector<SUPPORT *> Supports) const throw (MEDEXCEPTION)
2405 {
2406   const char * LOC = "MESH:::intersectSupports(const vector<SUPPORT *> ) : " ;
2407   BEGIN_OF(LOC) ;
2408
2409   SUPPORT * returnedSupport;
2410   string returnedSupportName;
2411   string returnedSupportDescription;
2412   char * returnedSupportNameChar;
2413   char * returnedSupportDescriptionChar;
2414   int size = Supports.size();
2415
2416   if (size == 1)
2417     {
2418       MESSAGE(LOC <<" there is only one SUPPORT in the argument list, the method return a copy of this object !");
2419       SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2420
2421       returnedSupport = new SUPPORT(*obj);
2422
2423       int lenName = strlen((Supports[0]->getName()).c_str()) + 8 + 1;
2424       int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 8 + 1;
2425
2426       returnedSupportNameChar = new char[lenName];
2427       returnedSupportDescriptionChar = new char[lenDescription];
2428
2429       returnedSupportNameChar = strcpy(returnedSupportNameChar,"Copy of ");
2430       returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
2431       returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Copy of ");
2432       returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2433                                               (Supports[0]->getDescription()).c_str());
2434
2435       returnedSupportName = string(returnedSupportNameChar);
2436       returnedSupportDescription = string(returnedSupportDescriptionChar);
2437
2438       returnedSupport->setName(returnedSupportName);
2439       returnedSupport->setDescription(returnedSupportDescription);
2440
2441       delete [] returnedSupportNameChar;
2442       delete [] returnedSupportDescriptionChar;
2443     }
2444   else
2445     {
2446       SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2447       returnedSupport = new SUPPORT(*obj);
2448
2449       int lenName = strlen((Supports[0]->getName()).c_str()) + 16 + 1;
2450       int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 16 + 1;
2451
2452       for(int i = 1;i<size;i++)
2453         {
2454           obj = const_cast <SUPPORT *> (Supports[i]);
2455           returnedSupport->intersecting(obj);
2456
2457           if (i == (size-1))
2458             {
2459               lenName = lenName + 5 + strlen((Supports[i]->getName()).c_str());
2460               lenDescription = lenDescription + 5 +
2461                 strlen((Supports[i]->getDescription()).c_str());
2462             }
2463           else
2464             {
2465               lenName = lenName + 2 + strlen((Supports[i]->getName()).c_str());
2466               lenDescription = lenDescription + 2 +
2467                 strlen((Supports[i]->getDescription()).c_str());
2468             }
2469         }
2470
2471       if(returnedSupport != (SUPPORT *) NULL)
2472         {
2473           returnedSupportNameChar = new char[lenName];
2474           returnedSupportDescriptionChar = new char[lenDescription];
2475
2476           returnedSupportNameChar = strcpy(returnedSupportNameChar,
2477                                            "Intersection of ");
2478           returnedSupportDescriptionChar =
2479             strcpy(returnedSupportDescriptionChar,"Intersection of ");
2480
2481           returnedSupportNameChar = strcat(returnedSupportNameChar,
2482                                            (Supports[0]->getName()).c_str());
2483           returnedSupportDescriptionChar =
2484             strcat(returnedSupportDescriptionChar,
2485                    (Supports[0]->getDescription()).c_str());
2486
2487           for(int i = 1;i<size;i++)
2488             {
2489               if (i == (size-1))
2490                 {
2491                   returnedSupportNameChar = strcat(returnedSupportNameChar,
2492                                                    " and ");
2493                   returnedSupportDescriptionChar =
2494                     strcat(returnedSupportDescriptionChar," and ");
2495
2496                   returnedSupportNameChar =
2497                     strcat(returnedSupportNameChar,
2498                            (Supports[i]->getName()).c_str());
2499                   returnedSupportDescriptionChar =
2500                     strcat(returnedSupportDescriptionChar,
2501                            (Supports[i]->getDescription()).c_str());
2502                 }
2503               else
2504                 {
2505                   returnedSupportNameChar = strcat(returnedSupportNameChar,
2506                                                    ", ");
2507                   returnedSupportNameChar =
2508                     strcat(returnedSupportNameChar,
2509                            (Supports[i]->getName()).c_str());
2510
2511                   returnedSupportDescriptionChar =
2512                     strcat(returnedSupportDescriptionChar,", ");
2513                   returnedSupportDescriptionChar =
2514                     strcat(returnedSupportDescriptionChar,
2515                            (Supports[i]->getDescription()).c_str());
2516                 }
2517             }
2518
2519           returnedSupportName = string(returnedSupportNameChar);
2520           returnedSupport->setName(returnedSupportName);
2521
2522           returnedSupportDescription = string(returnedSupportDescriptionChar);
2523           returnedSupport->setDescription(returnedSupportDescription);
2524
2525           delete [] returnedSupportNameChar;
2526           delete [] returnedSupportDescriptionChar;
2527         }
2528     }
2529
2530   END_OF(LOC) ;
2531   return returnedSupport;
2532 }
2533
2534
2535 // internal helper type
2536 struct _cell
2537 {
2538     //int number;
2539     std::vector<int> groups;
2540     MED_EN::medGeometryElement geometricType;
2541 };
2542
2543 // Create families from groups
2544 void MESH::createFamilies()
2545 {
2546     int nbFam=0; // count the families we create, used as identifier ???????????
2547
2548     int idFamNode = 0; // identifier for node families
2549     int idFamElement = 0; // identifier for cell, face or edge families
2550
2551     // Main loop on mesh's entities
2552     for (medEntityMesh entity=MED_CELL; entity!=MED_ALL_ENTITIES; ++entity)
2553     {
2554         int numberofgroups = getNumberOfGroups(entity);
2555         if(!numberofgroups)
2556             continue; // no groups for this entity
2557
2558         vector< vector<FAMILY*> > whichFamilyInGroup(numberofgroups); // this container is used to update groups at the end
2559
2560         // make myFamilies points to the member corresponding to entity
2561         vector<FAMILY*>* myFamilies;
2562         switch ( entity )
2563         {
2564             case MED_CELL :
2565                 myFamilies = & _familyCell;
2566                 break;
2567             case MED_FACE :
2568                 myFamilies = & _familyFace;
2569                 break;
2570             case MED_EDGE :
2571                 myFamilies = & _familyEdge;
2572                 break;
2573             case MED_NODE :
2574                 myFamilies = & _familyNode;
2575                 break;
2576         }
2577
2578         vector<GROUP*> myGroups=getGroups(entity); // get a copy of the groups ptr for the entity
2579         // get a copy of the (old) family ptrs before clearing
2580         vector<FAMILY*> myOldFamilies=getFamilies(entity);
2581         myFamilies->clear();
2582
2583
2584         // 1 - Create a vector containing for each cell (of the entity) an information structure
2585         //     giving geometric type and the groups it belong to
2586
2587         med_int numberOfTypes=getNumberOfTypes(entity);
2588         const int * index=getGlobalNumberingIndex(entity);
2589         const medGeometryElement* geometricTypes=_connectivity->getGeometricTypes(entity); // pb avec entity=MED_NODE???
2590         med_int numberOfCells=index[numberOfTypes]-1;  // total number of cells for that entity
2591         SCRUTE(numberOfTypes);
2592         SCRUTE(numberOfCells);
2593         vector< _cell > tab_cell(numberOfCells);
2594         for(med_int t=0; t!=numberOfTypes; ++t)
2595             for(int n=index[t]-1; n!=index[t+1]-1; ++n)
2596                 tab_cell[n].geometricType=geometricTypes[t];
2597
2598
2599         // 2 - Scan cells in groups and update in tab_cell the container of groups a cell belong to
2600
2601         for (unsigned g=0; g!=myGroups.size(); ++g)
2602         {
2603             // scan cells that belongs to the group
2604             const med_int* groupCells=myGroups[g]->getnumber()->getValue();
2605             int nbCells=myGroups[g]->getnumber()->getLength();
2606             for(int c=0; c!=nbCells; ++c)
2607                 tab_cell[groupCells[c]-1].groups.push_back(g);
2608         }
2609
2610
2611         // 3 - Scan the cells vector, genarate family name, and create a map associating the family names
2612         //     whith the vector of contained cells
2613
2614         map< string,vector<int> > tab_families;
2615         map< string,vector<int> >::iterator fam;
2616         for(int n=0; n!=numberOfCells; ++n)
2617         {
2618             ostringstream key; // to generate the name of the family
2619             key << "FAM";
2620             if(tab_cell[n].groups.empty()) // this cell don't belong to any group
2621                 key << "_NONE" << entity;
2622
2623             for(vector<int>::const_iterator it=tab_cell[n].groups.begin(); it!=tab_cell[n].groups.end(); ++it)
2624             {
2625                 string groupName=myGroups[*it]->getName();
2626                 if(groupName.empty())
2627                     key << "_G" << *it;
2628                 else
2629                     key << "_" << groupName;
2630             }
2631
2632             tab_families[key.str()].push_back(n+1); // fill the vector of contained cells associated whith the family
2633         /*    fam = tab_families.find(key.str());
2634             if( fam != tab_families.end())
2635                 fam->second.push_back(n+1); // +1 : convention Fortran de MED
2636             else
2637             {
2638                 vector<int> newfamily;
2639                 newfamily.push_back(n+1); // +1 : convention Fortran de MED
2640                 tab_families.insert(make_pair(key.str(),newfamily));
2641             }*/
2642         }
2643
2644
2645         // 4 - Scan the family map, create MED Families, check if it already exist.
2646
2647         for( fam=tab_families.begin(); fam!=tab_families.end(); ++fam)
2648         {
2649             vector<medGeometryElement> tab_types_geometriques;
2650             medGeometryElement geometrictype=MED_NONE;
2651             vector<int> tab_index_types_geometriques;
2652             vector<int> tab_nombres_elements;
2653
2654             // scan family cells and fill the tab that are needed by the create a MED FAMILY
2655             for( int i=0; i!=fam->second.size(); ++i)
2656             {
2657                 int ncell=fam->second[i]-1;
2658                 if(tab_cell[ncell].geometricType != geometrictype)
2659                 {
2660                     // new geometric type -> we store it and complete index tabs
2661                     if(!tab_index_types_geometriques.empty())
2662                         tab_nombres_elements.push_back(i+1-tab_index_types_geometriques.back());
2663                     tab_types_geometriques.push_back( (geometrictype=tab_cell[ncell].geometricType));
2664                     tab_index_types_geometriques.push_back(i+1);
2665                 }
2666             }
2667             // store and complete index tabs for the last geometric type
2668             tab_nombres_elements.push_back(fam->second.size()+1-tab_index_types_geometriques.back());
2669             tab_index_types_geometriques.push_back(fam->second.size()+1);
2670
2671             // create a empty MED FAMILY and fill it with the tabs we constructed
2672             FAMILY* newFam = new FAMILY();
2673             newFam->setTotalNumberOfElements(fam->second.size());
2674             newFam->setName(fam->first);
2675             newFam->setMesh(this);
2676             newFam->setNumberOfGeometricType(tab_types_geometriques.size());
2677             newFam->setGeometricType(&tab_types_geometriques[0]); // we know the tab is not empy
2678             newFam->setNumberOfElements(&tab_nombres_elements[0]);
2679             newFam->setNumber(&tab_index_types_geometriques[0],&fam->second[0]);
2680             newFam->setEntity(entity);
2681             newFam->setAll(false);
2682
2683             int idFam = 0;
2684
2685             switch ( entity )
2686               {
2687               case MED_NODE :
2688                 ++idFamNode;
2689                 idFam = idFamNode;
2690                 break;
2691               case MED_CELL:
2692                 ++idFamElement;
2693                 idFam = -idFamElement;
2694                 break;
2695               case MED_FACE :
2696                 ++idFamElement;
2697                 idFam = -idFamElement;
2698                 break;
2699               case MED_EDGE :
2700                 ++idFamElement;
2701                 idFam = -idFamElement;
2702                 break;
2703               }
2704
2705             newFam->setIdentifier(idFam);
2706
2707             // Update links between families and groups
2708
2709             int ncell1=fam->second[0]-1;  // number of first cell in family
2710             int numberOfGroups=tab_cell[ncell1].groups.size(); // number of groups in the family
2711             if(numberOfGroups)
2712             {
2713                 newFam->setNumberOfGroups(numberOfGroups);
2714                 string * groupNames=new string[numberOfGroups];
2715
2716                 // iterate on the groups the family belongs to
2717                 vector<int>::const_iterator it=tab_cell[ncell1].groups.begin();
2718                 for(int ng=0 ; it!=tab_cell[ncell1].groups.end(); ++it, ++ng)
2719                 {
2720                     whichFamilyInGroup[*it].push_back(newFam);
2721                     groupNames[ng]=myGroups[*it]->getName();
2722                 }
2723                 newFam->setGroupsNames(groupNames);
2724             }
2725
2726             int sizeOfFamVect = myFamilies->size();
2727
2728             MESSAGE("  MESH::createFamilies() entity " << entity << " size " << sizeOfFamVect);
2729
2730             myFamilies->push_back(newFam);
2731         }
2732
2733         // delete old families
2734         for (int i=0;i<myOldFamilies.size();i++)
2735             delete myOldFamilies[i] ;
2736
2737         // update references in groups
2738         for (int i=0;i<myGroups.size();i++)
2739         {
2740             myGroups[i]->setNumberOfFamilies(whichFamilyInGroup[i].size());
2741             myGroups[i]->setFamilies(whichFamilyInGroup[i]);
2742         }
2743
2744     // re-scan the cells vector, and fill the family vector with cells.
2745     // creation of support, check if it already exist.
2746     }
2747 }