src/MEDMEM/MEDMEM_Mesh.cxx

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