src/MEDMEM/MEDMEM_Mesh.cxx

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