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