src/MEDMEM/MEDMEM_Mesh.cxx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 /*
  24  File Mesh.cxx
  25 */
  26
  27 #include "MEDMEM_DriversDef.hxx"
  28 #include "MEDMEM_Field.hxx"
  29 #include "MEDMEM_Mesh.hxx"
  30
  31 #include "MEDMEM_Support.hxx"
  32 #include "MEDMEM_Family.hxx"
  33 #include "MEDMEM_Group.hxx"
  34 #include "MEDMEM_Coordinate.hxx"
  35 #include "MEDMEM_Connectivity.hxx"
  36 #include "MEDMEM_CellModel.hxx"
  37 #include "VolSurfFormulae.hxx"
  38 #include "MEDMEM_DriverFactory.hxx"
  39
  40 #include "PointLocator.hxx"
  41
  42 #include <math.h>
  43 #include <map>
  44 #include <sstream>
  45
  46 using namespace std;
  47 using namespace MEDMEM;
  48 using namespace MED_EN;
  49
  50 #define MED_NOPDT -1
  51 #define MED_NONOR -1
  52
  53 // Block defining groups for the MEDMEM_ug documentation
  54 /*!
  55 \defgroup MESH_constructors MESH Constructors
  56
  57 The MESH class provides only two constructors : a copy constructor and
  58 a constructor enabling creation from file reading. The creation of
  59 a user-defined mesh implies the use of the MESHING class.
  60
  61 \defgroup MESH_advanced MESH Advanced features
  62 These functions provide access to high-level manipulation of the meshes, giving
  63 information about the cells or extracting supports from the mesh.
  64
  65 \defgroup MESH_connectivity MESH Connectivity information
  66 These methods are related to the extraction of connectivity information
  67 from the mesh.
  68
  69 \defgroup MESH_nodes MESH Nodes information
  70 These methods are related to the extraction of information about the mesh nodes.
  71
  72 */
  73
  74 void MESH::init()
  75 {
  76   const char* LOC = "MESH::init(): ";
  77   BEGIN_OF_MED(LOC);
  78
  79   GMESH::init();
  80
  81   delete _coordinate;
  82   _coordinate   = (COORDINATE   *) NULL;
  83   delete _connectivity;
  84   _connectivity = (CONNECTIVITY *) NULL;
  85
  86   _numberOfNodes = MED_INVALID;
  87
  88   _arePresentOptionnalNodesNumbers = 0;
  89
  90   END_OF_MED(LOC);
  91 }
  92
  93 /*! Create an empty MESH. */
  94 MESH::MESH():GMESH(),_coordinate(NULL),_connectivity(NULL)
  95 {
  96   init();
  97 }
  98
  99 /*! \if MEDMEM_ug
 100   \addtogroup MESH_constructors
 101   @{
 102 \endif
 103 */
 104 /*!
 105   Copy constructor
 106 */
 107 MESH::MESH(MESH &m): GMESH(m)
 108 {
 109   if (m._coordinate != NULL)
 110     _coordinate = new COORDINATE(* m._coordinate);
 111   else
 112     _coordinate = (COORDINATE *) NULL;
 113
 114   if (m._connectivity != NULL)
 115     _connectivity = new CONNECTIVITY(* m._connectivity);
 116   else
 117     _connectivity = (CONNECTIVITY *) NULL;
 118
 119   _numberOfNodes = m._numberOfNodes;
 120
 121   _arePresentOptionnalNodesNumbers = m._arePresentOptionnalNodesNumbers;
 122   _optionnalToCanonicNodesNumbers = m._optionnalToCanonicNodesNumbers;
 123 }
 124
 125 /*!
 126 \if MEDMEM_ug
 127 @}
 128 \endif
 129 */
 130
 131 MESH::~MESH()
 132 {
 133   MESSAGE_MED("MESH::~MESH() : Destroying the Mesh");
 134
 135   if (_coordinate != ((COORDINATE *) NULL)) delete _coordinate ;
 136   if (_connectivity != ((CONNECTIVITY *) NULL)) delete _connectivity ;
 137   _coordinate = 0;
 138   _connectivity = 0;
 139 }
 140
 141 MESH & MESH::operator=(const MESH &m)
 142 {
 143   const char* LOC = "MESH & MESH::operator=(const MESH &m) : ";
 144   BEGIN_OF_MED(LOC);
 145
 146   _name = m._name;
 147   _description = m._description;
 148
 149   if ( _coordinate   ) delete _coordinate;
 150   _coordinate   = m._coordinate   ? new COORDINATE  ( *m._coordinate   ) : 0;
 151   if ( _connectivity ) delete _connectivity;
 152   _connectivity = m._connectivity ? new CONNECTIVITY( *m._connectivity ) : 0;
 153
 154   _spaceDimension = m._spaceDimension;
 155   _numberOfNodes  = m._numberOfNodes;
 156
 157   _arePresentOptionnalNodesNumbers = m._arePresentOptionnalNodesNumbers;
 158   _optionnalToCanonicNodesNumbers  = m._optionnalToCanonicNodesNumbers;
 159
 160   vector<FAMILY*>*        fams[4] = { &_familyNode, &_familyCell, &_familyFace, &_familyEdge};
 161   const vector<FAMILY*>* mfams[4] = { &m._familyNode,&m._familyCell,&m._familyFace,&m._familyEdge };
 162   for ( int i = 0; i < 4; ++i )
 163   {
 164     for ( int f = 0; f < fams[i]->size(); ++f )
 165       fams[i]->at(f)->removeReference();
 166     fams[i]->clear();
 167     fams[i]->reserve( mfams[i]->size() );
 168     for ( int f = 0; f < mfams[i]->size(); ++f )
 169     {
 170       if ( mfams[i]->at(f) )
 171       {
 172         fams[i]->push_back( new FAMILY( *mfams[i]->at(f) ));
 173         fams[i]->back()->setMesh( this );
 174       }
 175     }
 176   }
 177   vector<GROUP*>*        groups[4] = { &_groupNode, &_groupCell, &_groupFace, &_groupEdge };
 178   const vector<GROUP*>* mgroups[4] = { &m._groupNode, &m._groupCell, &m._groupFace, &m._groupEdge };
 179   for ( int i = 0; i < 4; ++i )
 180   {
 181     for ( int g = 0; g < groups[i]->size(); ++g )
 182       groups[i]->at(g)->removeReference();
 183     groups[i]->clear();
 184     groups[i]->reserve( mgroups[i]->size() );
 185     for ( int g = 0; g < mgroups[i]->size(); ++g )
 186     {
 187       if ( mgroups[i]->at(g) )
 188       {
 189         groups[i]->push_back( new GROUP( *mgroups[i]->at(g) ));
 190         groups[i]->back()->setMesh( this );
 191       }
 192     }
 193   }
 194
 195   for ( int drv = 0; drv < _drivers.size(); ++drv )
 196     delete _drivers[drv];
 197   _drivers.clear();
 198   _drivers.reserve( m._drivers.size());
 199   for ( int drv = 0; drv < m._drivers.size(); ++drv )
 200     if ( m._drivers[drv] )
 201       _drivers.push_back( m._drivers[drv]->copy() );
 202
 203   return *this;
 204 }
 205
 206 bool MESH::operator==(const MESH& other) const
 207 {
 208   const char* LOC = "MESH::operator==";
 209   BEGIN_OF_MED(LOC);
 210   return this==&other;
 211 }
 212
 213 /*!\if MEDMEM_ug
 214 \addtogroup MESH_constructors
 215 @{
 216 \endif
 217 */
 218 /*!
 219   Creates a %MESH object using a %MESH driver of type %driverTypes (MED_DRIVER, GIBI_DRIVER, ...) associated with file \a fileName. As several meshes can coexist in the same file (notably in MED files) , the constructor takes a third argument that specifies the name of the mesh.
 220   The constructor will throw an exception if the file does not exist, has no reading permissions or if the mesh does not exist in the file.
 221 */
 222 MESH::MESH(driverTypes driverType, const string &  fileName/*=""*/, const string &  driverName/*=""*/) throw (MEDEXCEPTION):_coordinate(0),_connectivity(0)
 223 {
 224   const char * LOC = "MESH::MESH(driverTypes driverType, const string &  fileName="", const string &  driverName/="") : ";
 225   BEGIN_OF_MED(LOC);
 226
 227   int current;
 228   init();
 229   GENDRIVER *myDriver=DRIVERFACTORY::buildDriverForMesh(driverType,fileName,this,driverName,RDONLY);
 230   current = addDriver(*myDriver);
 231   delete myDriver;
 232   try
 233   {
 234     _drivers[current]->open();
 235     _drivers[current]->read();
 236     _drivers[current]->close();
 237   }
 238   catch ( MED_EXCEPTION& e )
 239   {
 240     if ( _drivers[current] )
 241       _drivers[current]->close(), delete _drivers[current];
 242     _drivers[current] = 0;
 243     throw e;
 244   }
 245
 246   END_OF_MED(LOC);
 247 }
 248 /*!\if MEDMEM_ug
 249   @}
 250 \endif
 251 */
 252
 253 /*!
 254   Returns true if mesh \a other has same
 255   coordinates (to 1E-15 precision ) and same connectivity as the calling object.
 256   Information like name or description is not taken into account
 257   for the comparison.
 258 */
 259
 260 bool MESH::deepCompare(const GMESH& gother) const
 261 {
 262   if ( gother.getIsAGrid() != getIsAGrid())
 263     return false;
 264   const MESH& other = static_cast<const MESH&>( gother );
 265
 266   int size1=getSpaceDimension()*getNumberOfNodes();
 267   int size2=other.getSpaceDimension()*other.getNumberOfNodes();
 268   if(size1!=size2)
 269     return false;
 270
 271   const COORDINATE* CRD = other.getCoordinateptr();
 272   if( (!CRD && _coordinate) || (CRD && !(_coordinate)) ) {
 273     return false;
 274   }
 275
 276   bool ret=true;
 277   if( _coordinate ) {
 278     const double* coord1=getCoordinates(MED_FULL_INTERLACE);
 279     const double* coord2=other.getCoordinates(MED_FULL_INTERLACE);
 280     for(int i=0;i<size1 && ret;i++) {
 281       ret=(fabs(coord1[i]-coord2[i])<1e-15);
 282     }
 283   }
 284   if(ret) {
 285     const CONNECTIVITY* CNV = other.getConnectivityptr();
 286     if( (!CNV && _connectivity) || (CNV && !(_connectivity)) ) {
 287       return false;
 288     }
 289     if(_connectivity) {
 290       return _connectivity->deepCompare(*other._connectivity);
 291     }
 292   }
 293   return ret;
 294 }
 295
 296 /*!
 297  * \brief print my contents
 298  */
 299 ostream & ::MEDMEM::operator<<(ostream &os, const MESH &myMesh)
 300 {
 301   myMesh.printMySelf(os);
 302   return os;
 303 }
 304 void MESH::printMySelf(ostream &os) const
 305 {
 306   const MESH &myMesh = *this;
 307   int spacedimension = myMesh.getSpaceDimension();
 308   int meshdimension  = myMesh.getMeshDimension();
 309   int numberofnodes  = myMesh.getNumberOfNodes();
 310
 311   if ( spacedimension == MED_INVALID ) {
 312     os << " Empty mesh ...";
 313     return;
 314   }
 315
 316   os << "Space Dimension : " << spacedimension << endl << endl;
 317
 318   os << "Mesh Dimension : " << meshdimension << endl << endl;
 319
 320   if(myMesh.getCoordinateptr()) {
 321     const double * coordinates = myMesh.getCoordinates(MED_FULL_INTERLACE);
 322
 323     os << "SHOW NODES COORDINATES : " << endl;
 324     os << "Name :" << endl;
 325     const string * coordinatesnames = myMesh.getCoordinatesNames();
 326     if ( coordinatesnames ) {
 327       for(int i=0; i<spacedimension ; i++)
 328       {
 329         os << " - " << coordinatesnames[i] << endl;
 330       }
 331     }
 332     os << "Unit :" << endl;
 333     const string * coordinatesunits = myMesh.getCoordinatesUnits();
 334     if ( coordinatesunits ) {
 335       for(int i=0; i<spacedimension ; i++)
 336       {
 337         os << " - " << coordinatesunits[i] << endl;
 338       }
 339     }
 340     for(int i=0; i<numberofnodes ; i++)
 341     {
 342       os << "Nodes " << i+1 << " : ";
 343       for (int j=0; j<spacedimension ; j++)
 344         os << coordinates[i*spacedimension+j] << " ";
 345       os << endl;
 346     }
 347   }
 348
 349   if(myMesh.getConnectivityptr()) {
 350     os << endl << "SHOW CONNECTIVITY  :" << endl;
 351     os << *myMesh._connectivity << endl;
 352   }
 353
 354   medEntityMesh entity;
 355   os << endl << "SHOW FAMILIES :" << endl << endl;
 356   for (int k=1; k<=4; k++)
 357   {
 358     if (k==1) entity = MED_NODE;
 359     if (k==2) entity = MED_CELL;
 360     if (k==3) entity = MED_FACE;
 361     if (k==4) entity = MED_EDGE;
 362     int numberoffamilies = myMesh.getNumberOfFamilies(entity);
 363     os << "NumberOfFamilies on "<< entNames[entity] <<" : "<<numberoffamilies<<endl;
 364     using namespace MED_EN;
 365     for (int i=1; i<numberoffamilies+1;i++)
 366     {
 367       os << * myMesh.getFamily(entity,i) << endl;
 368     }
 369   }
 370
 371   os << endl << "SHOW GROUPS :" << endl << endl;
 372   for (int k=1; k<=4; k++)
 373   {
 374     if (k==1) entity = MED_NODE;
 375     if (k==2) entity = MED_CELL;
 376     if (k==3) entity = MED_FACE;
 377     if (k==4) entity = MED_EDGE;
 378     int numberofgroups = myMesh.getNumberOfGroups(entity);
 379     os << "NumberOfGroups on "<< entNames[entity] <<" : "<<numberofgroups<<endl;
 380     using namespace MED_EN;
 381     for (int i=1; i<numberofgroups+1;i++)
 382     {
 383       os << * myMesh.getGroup(entity,i) << endl;
 384     }
 385   }
 386 }
 387
 388 /*! \if MEDMEM_ug
 389 \addtogroup MESH_general
 390 @{
 391 \endif
 392 */
 393 /*! Gets the dimension of the mesh (2 for 2D- and 3D-surfaces, 3 for volumes). */
 394 int MESH::getMeshDimension() const
 395 {
 396   int dim = -1;
 397   if ( _connectivity )
 398     for ( int i = 0; i < _connectivity->getNumberOfTypes(MED_EN::MED_CELL); ++i )
 399       if ( _connectivity->getCellsTypes(MED_EN::MED_CELL)[i].getDimension() > dim )
 400         dim = _connectivity->getCellsTypes(MED_EN::MED_CELL)[i].getDimension();
 401   return dim;
 402 }
 403 /*! \if MEDMEM_ug @} \endif */
 404
 405 /*!
 406   Get global number of element which have same connectivity than connectivity argument.
 407
 408   It do not take care of connectivity order (3,4,7,10 is same as 7,3,10,4).
 409
 410   Return -1 if not found.
 411 */
 412 int MESH::getElementNumber(MED_EN::medConnectivity ConnectivityType,
 413                            MED_EN::medEntityMesh Entity,
 414                            MED_EN::medGeometryElement Type,
 415                            int * connectivity) const
 416 {
 417   const char* LOC="MESH::getElementNumber " ;
 418   BEGIN_OF_MED(LOC) ;
 419
 420   int numberOfValue ; // size of connectivity array
 421   CELLMODEL myType(Type) ;
 422   if (ConnectivityType==MED_DESCENDING)
 423     numberOfValue = myType.getNumberOfConstituents(1) ; // faces (3D) or edges (2D)
 424   else
 425     numberOfValue = myType.getNumberOfNodes() ; // nodes
 426
 427   const int * myReverseConnectivityValue = getReverseConnectivity(ConnectivityType,Entity) ;
 428   const int * myReverseConnectivityIndex = getReverseConnectivityIndex(ConnectivityType,Entity) ;
 429
 430   // First node or face/edge
 431   int indexBegin = myReverseConnectivityIndex[connectivity[0]-1] ;
 432   int indexEnd = myReverseConnectivityIndex[connectivity[0]] ;
 433
 434   // list to put cells numbers
 435   list<int> cellsList ;
 436   list<int>::iterator itList ;
 437   for (int i=indexBegin; i<indexEnd; i++)
 438     cellsList.push_back(myReverseConnectivityValue[i-1]) ;
 439
 440   // Foreach node or face/edge in cell, we search cells which are in common.
 441   // At the end, we must have only one !
 442
 443   for (int i=1; i<numberOfValue; i++) {
 444     int connectivity_i = connectivity[i] ;
 445     for (itList=cellsList.begin();itList!=cellsList.end();/*itList++*/) {
 446       bool find = false ;
 447       for (int j=myReverseConnectivityIndex[connectivity_i-1]; j<myReverseConnectivityIndex[connectivity_i]; j++) {
 448         if ((*itList)==myReverseConnectivityValue[j-1]) {
 449           find=true ;
 450           break ;
 451         }
 452       }
 453       if (!find)
 454         itList=cellsList.erase(itList++);
 455       else
 456         itList++;
 457     }
 458   }
 459
 460   if (cellsList.size()>1) // we have more than one cell
 461     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Found more than one element !")) ;
 462
 463   if (cellsList.size()==0)
 464     return -1;
 465
 466   END_OF_MED(LOC);
 467
 468   return cellsList.front() ;
 469 }
 470
 471 /*!
 472 \addtogroup MESH_advanced
 473 @{
 474 The methods described in this section are algorithms that perform a computation
 475 and return a result in the form of a SUPPORT or a FIELD to the user. For large meshes,
 476 as the returned information is not stored in the mesh but is rather computed, the
 477 computation time can be large.
 478 */
 479
 480 /*!
 481   Returns a support which reference all elements on the boundary of mesh.
 482   For a d-dimensional mesh, a boundary element is defined as a d-1 dimension
 483   element that is referenced by only one element in the full descending connectivity.
 484
 485   This method can also return the list of nodes that belong to the boundary elements.
 486
 487   WARNING: This method can recalculate descending connectivity from partial to full form,
 488   so that partial SUPPORT on d-1 dimension elements becomes invalid.
 489
 490   \param Entity entity on which the boundary is desired. It has to be either \a MED_NODE or the
 491   d-1 dimension entity type (MED_FACE in 3D, MED_EDGE in 2D).
 492 */
 493 SUPPORT * MESH::getBoundaryElements(MED_EN::medEntityMesh Entity) const throw (MEDEXCEPTION)
 494 {
 495   const char * LOC = "MESH::getBoundaryElements : " ;
 496   BEGIN_OF_MED(LOC) ;
 497   // some test :
 498   // actually we could only get face (in 3D) and edge (in 2D)
 499   medEntityMesh entityToParse=Entity;
 500   if(_spaceDimension == 3)
 501     if (Entity != MED_FACE)
 502       {
 503         if(Entity==MED_NODE)
 504           entityToParse=MED_FACE;
 505         else
 506           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<Entity<<" !"));
 507       }
 508   if(_spaceDimension == 2)
 509     if(Entity != MED_EDGE)
 510       {
 511         if(Entity==MED_NODE)
 512           entityToParse=MED_EDGE;
 513         else
 514           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<Entity<<" !"));
 515       }
 516
 517   // assure that descending connectivity is full
 518   if ( !_connectivity )
 519     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<< "no connectivity defined in MESH !"));
 520   _connectivity->calculateFullDescendingConnectivity(MED_CELL);
 521
 522   const int * myConnectivityValue = getReverseConnectivity(MED_DESCENDING) ;
 523   const int * myConnectivityIndex = getReverseConnectivityIndex(MED_DESCENDING) ;
 524   int numberOf = getNumberOfElements(entityToParse,MED_ALL_ELEMENTS) ;
 525   list<int> myElementsList;
 526
 527   for (int i=0 ; i<numberOf; i++)
 528     if (myConnectivityValue[myConnectivityIndex[i]] == 0)
 529       myElementsList.push_back(i+1);
 530
 531   if ( myElementsList.empty() && numberOf != 0 )
 532     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No boundary elements found by reverse descending connectivity for entity "<<Entity<<" !"));
 533
 534   if(Entity==MED_NODE)
 535     return buildSupportOnNodeFromElementList(myElementsList,entityToParse);
 536   else
 537     return buildSupportOnElementsFromElementList(myElementsList,entityToParse);
 538 }
 539 /*!
 540   @}
 541 */
 542
 543 /*!
 544   Method that do the same thing as buildSupportOnNodeFromElementList except that a SUPPORT is not created.
 545 */
 546 void MESH::fillSupportOnNodeFromElementList(const list<int>& listOfElt, SUPPORT *supportToFill) const throw (MEDEXCEPTION)
 547 {
 548   MED_EN::medEntityMesh entity=supportToFill->getEntity();
 549   supportToFill->setAll(false);
 550   supportToFill->setMesh((MESH *)this);
 551
 552   set<int> nodes;
 553   if ( entity == MED_NODE )
 554     {
 555       supportToFill->fillFromNodeList(listOfElt);
 556     }
 557   else
 558     {
 559       int lgth;
 560       for(list<int>::const_iterator iter=listOfElt.begin();iter!=listOfElt.end();iter++)
 561         {
 562           const int *conn=_connectivity->getConnectivityOfAnElement(MED_NODAL,entity,*iter,lgth);
 563           nodes.insert( conn, conn+lgth );
 564         }
 565       list<int> nodesList( nodes.begin(), nodes.end() );
 566       supportToFill->fillFromNodeList( nodesList );
 567     }
 568 }
 569
 570 /*!
 571   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
 572   entity 'entity'.
 573 */
 574 SUPPORT *MESH::buildSupportOnNodeFromElementList(const list<int>& listOfElt,MED_EN::medEntityMesh entity) const throw (MEDEXCEPTION)
 575 {
 576   SUPPORT * mySupport = new SUPPORT;
 577   mySupport->setMesh((MESH *)this);
 578   mySupport->setName("Boundary");
 579   mySupport->setEntity( entity );
 580   fillSupportOnNodeFromElementList(listOfElt,mySupport);
 581   return mySupport;
 582 }
 583
 584 /*!
 585 \addtogroup MESH_advanced
 586 @{
 587 */
 588 /*! Retrieves the volume of all the elements contained in \a Support. This method returns
 589   a FIELD structure based on this support. It only works on MED_CELL for 3D meshes.
 590 */
 591 FIELD<double, FullInterlace>* MESH::getVolume(const SUPPORT *Support, bool isAbs) const throw (MEDEXCEPTION)
 592 {
 593   const char * LOC = "MESH::getVolume(SUPPORT*) : ";
 594   BEGIN_OF_MED(LOC);
 595
 596   // Support must be on 3D elements
 597
 598   // Make sure that the MESH class is the same as the MESH class attribut
 599   // in the class Support
 600   const GMESH* myMesh = Support->getMesh();
 601   if (this != myMesh)
 602     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
 603
 604   int dim_space = getSpaceDimension();
 605   medEntityMesh support_entity = Support->getEntity();
 606   bool onAll = Support->isOnAllElements();
 607
 608   int nb_type, length_values;
 609   const medGeometryElement* types;
 610   int nb_entity_type;
 611   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
 612   const int* global_connectivity;
 613   nb_type = Support->getNumberOfTypes();
 614   length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
 615   types = Support->getTypes();
 616   int index;
 617   FIELD<double, FullInterlace>* Volume =
 618     new FIELD<double, FullInterlace>(Support,1);
 619   Volume->setName("VOLUME");
 620   Volume->setDescription("cells volume");
 621   Volume->setComponentName(1,"volume");
 622   Volume->setComponentDescription(1,"desc-comp");
 623
 624   string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1]+"*"+myMesh->getCoordinatesUnits()[2];
 625
 626   Volume->setMEDComponentUnit(1,MEDComponentUnit);
 627
 628   Volume->setIterationNumber(0);
 629   Volume->setOrderNumber(0);
 630   Volume->setTime(0.0);
 631
 632   typedef  MEDMEM_ArrayInterface<double,FullInterlace,NoGauss>::Array ArrayNoGauss;
 633   ArrayNoGauss  *volume = Volume->getArrayNoGauss();
 634
 635   index = 1;
 636   const double * coord = getCoordinates(MED_FULL_INTERLACE);
 637
 638   for (int i=0;i<nb_type;i++)
 639   {
 640     medGeometryElement type = types[i] ;
 641     double xvolume;
 642     nb_entity_type = Support->getNumberOfElements(type);
 643     if(type != MED_EN::MED_POLYHEDRA)
 644     {
 645       if (onAll)
 646       {
 647         global_connectivity = getConnectivity(MED_NODAL,support_entity,type);
 648       }
 649       else
 650       {
 651         const int * supp_number = Support->getNumber(type);
 652         const int * connectivity = getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
 653         const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
 654         int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
 655
 656         for (int k_type = 0; k_type<nb_entity_type; k_type++) {
 657           for (int j_ent = 0; j_ent<(type%100); j_ent++) {
 658             global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
 659           }
 660         }
 661         global_connectivity = global_connectivity_tmp ;
 662       }
 663     }
 664
 665     switch (type)
 666     {
 667     case MED_TETRA4 : case MED_TETRA10 :
 668       {
 669         for (int tetra=0;tetra<nb_entity_type;tetra++)
 670         {
 671           int tetra_index = (type%100)*tetra;
 672           int N1 = global_connectivity[tetra_index]-1;
 673           int N2 = global_connectivity[tetra_index+1]-1;
 674           int N3 = global_connectivity[tetra_index+2]-1;
 675           int N4 = global_connectivity[tetra_index+3]-1;
 676           xvolume=INTERP_KERNEL::calculateVolumeForTetra(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4);
 677           if(isAbs)
 678             xvolume=fabs(xvolume);
 679           volume->setIJ(index,1,xvolume) ;
 680           index++;
 681         }
 682         break;
 683       }
 684     case MED_PYRA5 : case MED_PYRA13 :
 685       {
 686         for (int pyra=0;pyra<nb_entity_type;pyra++)
 687         {
 688           int pyra_index = (type%100)*pyra;
 689           int N1 = global_connectivity[pyra_index]-1;
 690           int N2 = global_connectivity[pyra_index+1]-1;
 691           int N3 = global_connectivity[pyra_index+2]-1;
 692           int N4 = global_connectivity[pyra_index+3]-1;
 693           int N5 = global_connectivity[pyra_index+4]-1;
 694           xvolume=INTERP_KERNEL::calculateVolumeForPyra(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4,coord+dim_space*N5);
 695           if(isAbs)
 696             xvolume=fabs(xvolume);
 697           volume->setIJ(index,1,xvolume) ;
 698           index++;
 699         }
 700         break;
 701       }
 702     case MED_PENTA6 : case MED_PENTA15 :
 703       {
 704         for (int penta=0;penta<nb_entity_type;penta++)
 705         {
 706           int penta_index = (type%100)*penta;
 707           int N1 = global_connectivity[penta_index]-1;
 708           int N2 = global_connectivity[penta_index+1]-1;
 709           int N3 = global_connectivity[penta_index+2]-1;
 710           int N4 = global_connectivity[penta_index+3]-1;
 711           int N5 = global_connectivity[penta_index+4]-1;
 712           int N6 = global_connectivity[penta_index+5]-1;
 713           xvolume=INTERP_KERNEL::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);
 714           if(isAbs)
 715             xvolume=fabs(xvolume);
 716           volume->setIJ(index,1,xvolume) ;
 717           index++;
 718         }
 719         break;
 720       }
 721     case MED_HEXA8 : case MED_HEXA20 :
 722       {
 723         for (int hexa=0;hexa<nb_entity_type;hexa++)
 724         {
 725           int hexa_index = (type%100)*hexa;
 726
 727           int N1 = global_connectivity[hexa_index]-1;
 728           int N2 = global_connectivity[hexa_index+1]-1;
 729           int N3 = global_connectivity[hexa_index+2]-1;
 730           int N4 = global_connectivity[hexa_index+3]-1;
 731           int N5 = global_connectivity[hexa_index+4]-1;
 732           int N6 = global_connectivity[hexa_index+5]-1;
 733           int N7 = global_connectivity[hexa_index+6]-1;
 734           int N8 = global_connectivity[hexa_index+7]-1;
 735           xvolume=INTERP_KERNEL::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);
 736           if(isAbs)
 737             xvolume=fabs(xvolume);
 738           volume->setIJ(index,1,xvolume) ;
 739           index++;
 740         }
 741         break;
 742       }
 743     case MED_POLYHEDRA:
 744       {
 745         double bary[3];
 746         if(onAll)
 747         {
 748           for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
 749           {
 750             int lgthNodes,iPts,iFaces,iPtsInFace;
 751             int offsetWithClassicType=getGlobalNumberingIndex(support_entity)[i]-1;
 752             int *nodes=_connectivity->getNodesOfPolyhedron(offsetWithClassicType+polyhs+1,lgthNodes);
 753             int nbOfFaces,*nbOfNodesPerFaces;
 754             int **nodes1=_connectivity->getNodesPerFaceOfPolyhedron(offsetWithClassicType+polyhs+1,nbOfFaces,nbOfNodesPerFaces);
 755             double **pts=new double * [lgthNodes];
 756             double ***pts1=new double ** [nbOfFaces];
 757             for(iPts=0;iPts<lgthNodes;iPts++)
 758               pts[iPts]=(double *)(coord+3*(nodes[iPts]-1));
 759             for(iFaces=0;iFaces<nbOfFaces;iFaces++)
 760             {
 761               pts1[iFaces]=new double* [nbOfNodesPerFaces[iFaces]];
 762               for(iPtsInFace=0;iPtsInFace<nbOfNodesPerFaces[iFaces];iPtsInFace++)
 763                 pts1[iFaces][iPtsInFace]=(double *)(coord+3*(nodes1[iFaces][iPtsInFace]-1));
 764             }
 765             delete [] nodes1;
 766             INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,lgthNodes,3,bary);
 767             delete [] nodes;
 768             delete [] pts;
 769             if(isAbs)
 770               xvolume=INTERP_KERNEL::calculateVolumeForPolyhAbs((const double ***)pts1,nbOfNodesPerFaces,nbOfFaces,bary);
 771             else
 772               xvolume=INTERP_KERNEL::calculateVolumeForPolyh((const double ***)pts1,nbOfNodesPerFaces,nbOfFaces,bary);
 773             delete [] nbOfNodesPerFaces;
 774             for(iFaces=0;iFaces<nbOfFaces;iFaces++)
 775               delete [] pts1[iFaces];
 776             delete [] pts1;
 777             volume->setIJ(index,1,xvolume) ;
 778             index++;
 779           }
 780         }
 781         else
 782         {
 783           const int * supp_number = Support->getNumber(MED_EN::MED_POLYHEDRA);
 784           for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
 785           {
 786             int lgthNodes,iPts,iFaces,iPtsInFace;
 787             int *nodes=_connectivity->getNodesOfPolyhedron(supp_number[polyhs],lgthNodes);
 788             int nbOfFaces,*nbOfNodesPerFaces;
 789             int **nodes1=_connectivity->getNodesPerFaceOfPolyhedron(supp_number[polyhs],nbOfFaces,nbOfNodesPerFaces);
 790             double **pts=new double * [lgthNodes];
 791             double ***pts1=new double ** [nbOfFaces];
 792             for(iPts=0;iPts<lgthNodes;iPts++)
 793               pts[iPts]=(double *)(coord+3*(nodes[iPts]-1));
 794             for(iFaces=0;iFaces<nbOfFaces;iFaces++)
 795             {
 796               pts1[iFaces]=new double* [nbOfNodesPerFaces[iFaces]];
 797               for(iPtsInFace=0;iPtsInFace<nbOfNodesPerFaces[iFaces];iPtsInFace++)
 798                 pts1[iFaces][iPtsInFace]=(double *)(coord+3*(nodes1[iFaces][iPtsInFace]-1));
 799             }
 800             delete [] nodes1;
 801             INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,lgthNodes,3,bary);
 802             delete [] nodes;
 803             delete [] pts;
 804             if(isAbs)
 805               xvolume=INTERP_KERNEL::calculateVolumeForPolyhAbs((const double ***)pts1,nbOfNodesPerFaces,nbOfFaces,bary);
 806             else
 807               xvolume=INTERP_KERNEL::calculateVolumeForPolyh((const double ***)pts1,nbOfNodesPerFaces,nbOfFaces,bary);
 808             delete [] nbOfNodesPerFaces;
 809             for(iFaces=0;iFaces<nbOfFaces;iFaces++)
 810               delete [] pts1[iFaces];
 811             delete [] pts1;
 812             volume->setIJ(index,1,xvolume) ;
 813             index++;
 814           }
 815         }
 816         break;
 817       }
 818     default :
 819       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Volumes on it !"));
 820       break;
 821     }
 822
 823     if (!onAll && type!=MED_EN::MED_POLYHEDRA)
 824       delete [] global_connectivity ;
 825   }
 826
 827   return Volume;
 828 }
 829 /*! Retrieves the area of all the elements contained in \a Support. This method returns
 830   a FIELD structure based on this support. It only works on MED_CELL for 2D meshes or MED_FACE
 831   for 3D meshes.
 832 */
 833
 834 FIELD<double, FullInterlace>* MESH::getArea(const SUPPORT * Support) const throw (MEDEXCEPTION)
 835 {
 836   const char * LOC = "MESH::getArea(SUPPORT*) : ";
 837   BEGIN_OF_MED(LOC);
 838
 839   // Support must be on 2D elements
 840
 841   // Make sure that the MESH class is the same as the MESH class attribut
 842   // in the class Support
 843   const GMESH* myMesh = Support->getMesh();
 844   if (this != myMesh)
 845     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
 846
 847   int dim_space = getSpaceDimension();
 848   medEntityMesh support_entity = Support->getEntity();
 849   bool onAll = Support->isOnAllElements();
 850
 851   int nb_type, length_values;
 852   const medGeometryElement* types;
 853   int nb_entity_type;
 854   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
 855   const int* global_connectivity;
 856
 857   nb_type = Support->getNumberOfTypes();
 858   length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
 859   types = Support->getTypes();
 860
 861   int index;
 862   FIELD<double, FullInterlace>* Area;
 863
 864   Area = new FIELD<double, FullInterlace>(Support,1);
 865   Area->setName("AREA");
 866   Area->setDescription("cells or faces area");
 867
 868   Area->setComponentName(1,"area");
 869   Area->setComponentDescription(1,"desc-comp");
 870
 871   /*  string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1];*/
 872
 873   string MEDComponentUnit="trucmuch";
 874
 875   Area->setMEDComponentUnit(1,MEDComponentUnit);
 876
 877   Area->setIterationNumber(0);
 878   Area->setOrderNumber(0);
 879   Area->setTime(0.0);
 880
 881   double *area = (double *)Area->getValue();
 882
 883   const double * coord = getCoordinates(MED_FULL_INTERLACE);
 884   index = 0;
 885
 886   for (int i=0;i<nb_type;i++)
 887   {
 888     medGeometryElement type = types[i] ;
 889     nb_entity_type = Support->getNumberOfElements(type);
 890     const int *global_connectivityIndex = 0;
 891     if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
 892     {
 893       global_connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
 894       if (onAll)
 895       {
 896         global_connectivity = getConnectivity(MED_NODAL,support_entity,type);
 897       }
 898       else
 899       {
 900         const int * supp_number = Support->getNumber(type);
 901         const int * connectivity = getConnectivity(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]=INTERP_KERNEL::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]=INTERP_KERNEL::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           int offsetWithClassicType=getGlobalNumberingIndex(support_entity)[i]-1;
 955           const int * connectivity = _connectivity->getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
 956           const int * connectivity_index = _connectivity->getConnectivityIndex(MED_EN::MED_NODAL,support_entity) + offsetWithClassicType;
 957           for (int polygs=0;polygs<nb_entity_type;polygs++)
 958           {
 959             int size=connectivity_index[polygs+1]-connectivity_index[polygs];
 960             double **pts=new double * [size];
 961             for(int iPts=0;iPts<size;iPts++)
 962               pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[polygs]+iPts-1]-1));
 963             area[index] = INTERP_KERNEL::calculateAreaForPolyg((const double **)pts,size,dim_space);
 964             delete [] pts;
 965             index++;
 966           }
 967         }
 968         else
 969         {
 970           const int * supp_number = Support->getNumber(MED_EN::MED_POLYGON);
 971           const int * connectivity = _connectivity->getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
 972           const int * connectivity_index = _connectivity->getConnectivityIndex(MED_EN::MED_NODAL,support_entity);
 973           for (int polygs=0;polygs<nb_entity_type;polygs++)
 974           {
 975             int size=connectivity_index[supp_number[polygs]]-connectivity_index[supp_number[polygs]-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]-1]+iPts-1]-1));
 979             area[index]=INTERP_KERNEL::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 /*! Retrieves the length of all the elements contained in \a Support. This method returns
 998   a FIELD structure based on this support. It only works on MED_EDGE supports.
 999 */
1000 FIELD<double, FullInterlace>* MESH::getLength(const SUPPORT * Support) const throw (MEDEXCEPTION)
1001 {
1002   const char * LOC = "MESH::getLength(SUPPORT*) : ";
1003   BEGIN_OF_MED(LOC);
1004
1005   // Support must be on 1D elements
1006
1007   // Make sure that the MESH class is the same as the MESH class attribut
1008   // in the class Support
1009   const GMESH* myMesh = Support->getMesh();
1010   if (this != myMesh)
1011     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1012
1013   int dim_space = getSpaceDimension();
1014   medEntityMesh support_entity = Support->getEntity();
1015   bool onAll = Support->isOnAllElements();
1016
1017   int nb_type, length_values;
1018   const medGeometryElement* types;
1019   int nb_entity_type;
1020   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1021   const int* global_connectivity;
1022
1023   nb_type = Support->getNumberOfTypes();
1024   length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1025   types = Support->getTypes();
1026
1027   int index;
1028   FIELD<double, FullInterlace>* Length;
1029
1030   Length = new FIELD<double, FullInterlace>(Support,1);
1031
1032   typedef  MEDMEM_ArrayInterface<double,FullInterlace,NoGauss>::Array ArrayNoGauss;
1033   ArrayNoGauss * length = Length->getArrayNoGauss();
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_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_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 /*! Retrieves the normal for all elements contained in SUPPORT \a Support.
1104   The method is only functional for 2D supports for 3D meshes and 1D supports
1105   for 2D meshes. It returns
1106   a FIELD for which the number of components is equal to the dimension of the
1107   mesh and which represents coordinates of the vector normal to the element.
1108   The direction of the vector is undetermined.
1109 */
1110 FIELD<double, FullInterlace>* MESH::getNormal(const SUPPORT * Support) const throw (MEDEXCEPTION)
1111 {
1112   const char * LOC = "MESH::getNormal(SUPPORT*) : ";
1113   BEGIN_OF_MED(LOC);
1114
1115   // Support must be on 2D or 1D elements
1116
1117   // Make sure that the MESH class is the same as the MESH class attribut
1118   // in the class Support
1119   const GMESH* myMesh = Support->getMesh();
1120   if (this != myMesh)
1121     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : pointeur problem !"));
1122
1123   int dim_space = getSpaceDimension();
1124   int mesh_dim=getMeshDimension();
1125   medEntityMesh support_entity = Support->getEntity();
1126   bool onAll = Support->isOnAllElements();
1127
1128   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 ))
1129     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"incompatible mesh dimension and entity"));
1130   int nb_type, length_values;
1131   const medGeometryElement* types;
1132   int nb_entity_type;
1133   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1134   const int* global_connectivity;
1135
1136   nb_type = Support->getNumberOfTypes();
1137   length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1138   types = Support->getTypes();
1139
1140   int index;
1141
1142   FIELD<double, FullInterlace>* Normal =
1143     new FIELD<double, FullInterlace>(Support,dim_space);
1144   Normal->setName("NORMAL");
1145   Normal->setDescription("cells or faces normal");
1146   for (int k=1;k<=dim_space;k++) {
1147     Normal->setComponentName(k,"normal");
1148     Normal->setComponentDescription(k,"desc-comp");
1149     Normal->setMEDComponentUnit(k,"unit");
1150   }
1151
1152   Normal->setIterationNumber(MED_NOPDT);
1153   Normal->setOrderNumber(MED_NONOR);
1154   Normal->setTime(0.0);
1155   double *normal = (double *)Normal->getValue();
1156
1157   const double * coord = getCoordinates(MED_FULL_INTERLACE);
1158   index = 0;
1159
1160   for (int i=0;i<nb_type;i++)
1161   {
1162     medGeometryElement type = types[i] ;
1163     nb_entity_type = Support->getNumberOfElements(type);
1164
1165     // Make sure we trying to get Normals on
1166     // TRIA3 or TRIA6 or QUAD4 QUAD8 (2D) cells on a 3D mesh
1167     // or on SEG2 or SEG3 (1D) cells on a 2D mesh
1168
1169     if ( (((type==MED_TRIA3) || (type==MED_TRIA6) ||
1170            (type==MED_QUAD4) || (type==MED_QUAD8) || (type==MED_POLYGON)) &&
1171           (dim_space != 3)) || (((type==MED_SEG2) || (type==MED_SEG3)) &&
1172                                 (dim_space != 2)) )
1173       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : dimension problem !"));
1174     if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
1175     {
1176       if (onAll)
1177       {
1178         global_connectivity = getConnectivity(MED_NODAL,support_entity,type);
1179       }
1180       else
1181       {
1182         const int * supp_number = Support->getNumber(type);
1183         const int * connectivity = getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1184         const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1185         int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1186
1187         for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1188           for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1189             global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1190           }
1191         }
1192
1193         global_connectivity = global_connectivity_tmp ;
1194       }
1195     }
1196
1197     switch (type)
1198     {
1199     case MED_TRIA3 : case MED_TRIA6 :
1200       {
1201         for (int tria=0;tria<nb_entity_type;tria++)
1202         {
1203           int tria_index = (type%100)*tria;
1204           int N1 = global_connectivity[tria_index]-1;
1205           int N2 = global_connectivity[tria_index+1]-1;
1206           int N3 = global_connectivity[tria_index+2]-1;
1207           INTERP_KERNEL::calculateNormalForTria(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,normal+3*index);
1208           index++;
1209         }
1210         break;
1211       }
1212     case MED_QUAD4 : case MED_QUAD8 :
1213       {
1214         for (int quad=0;quad<nb_entity_type;quad++)
1215         {
1216           int quad_index = (type%100)*quad;
1217           int N1 = global_connectivity[quad_index]-1;
1218           int N2 = global_connectivity[quad_index+1]-1;
1219           int N3 = global_connectivity[quad_index+2]-1;
1220           int N4 = global_connectivity[quad_index+3]-1;
1221           INTERP_KERNEL::calculateNormalForQuad(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4,normal+3*index);
1222           index++;
1223         }
1224         break;
1225       }
1226     case MED_SEG2 : case MED_SEG3 :
1227       {
1228         double xnormal1, xnormal2;
1229         for (int edge=0;edge<nb_entity_type;edge++)
1230         {
1231           int edge_index = (type%100)*edge;
1232
1233           int N1 = global_connectivity[edge_index]-1;
1234           int N2 = global_connectivity[edge_index+1]-1;
1235
1236           double x1 = coord[dim_space*N1];
1237           double x2 = coord[dim_space*N2];
1238
1239           double y1 = coord[dim_space*N1+1];
1240           double y2 = coord[dim_space*N2+1];
1241
1242           xnormal1 = -(y2-y1);
1243           xnormal2 = x2-x1;
1244
1245           normal[2*index] = xnormal1 ;
1246           normal[2*index+1] = xnormal2 ;
1247
1248           index++;
1249         }
1250         break;
1251       }
1252     case MED_POLYGON :
1253       {
1254         if(onAll)
1255         {
1256           int offsetWithClassicType=getGlobalNumberingIndex(support_entity)[i]-1;
1257           const int * connectivity = _connectivity->getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1258           const int * connectivity_index = _connectivity->getConnectivityIndex(MED_EN::MED_NODAL,support_entity) + offsetWithClassicType;
1259           for (int polygs=0;polygs<nb_entity_type;polygs++)
1260           {
1261             int size=connectivity_index[polygs+1]-connectivity_index[polygs];
1262             double **pts=new double * [size];
1263             for(int iPts=0;iPts<size;iPts++)
1264               pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[polygs]+iPts-1])-1);
1265             INTERP_KERNEL::calculateNormalForPolyg((const double **)pts,size,normal+3*index);
1266             delete [] pts;
1267             index++;
1268           }
1269         }
1270         else
1271         {
1272           const int * supp_number = Support->getNumber(MED_EN::MED_POLYGON);
1273           const int * connectivity = _connectivity->getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1274           const int * connectivity_index = _connectivity->getConnectivityIndex(MED_EN::MED_NODAL,support_entity);
1275           for (int polygs=0;polygs<nb_entity_type;polygs++)
1276           {
1277             int size=connectivity_index[supp_number[polygs]]-connectivity_index[supp_number[polygs]-1];
1278             double **pts=new double * [size];
1279             for(int iPts=0;iPts<size;iPts++)
1280               pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[supp_number[polygs]-1]+iPts-1])-1);
1281             INTERP_KERNEL::calculateNormalForPolyg((const double **)pts,size,normal+3*index);
1282             delete [] pts;
1283             index++;
1284           }
1285         }
1286         break;
1287       }
1288     default :
1289       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Normals on it !"));
1290       break;
1291     }
1292     if (!onAll && type!=MED_EN::MED_POLYGON)
1293       delete [] global_connectivity ;
1294   }
1295   END_OF_MED(LOC);
1296
1297   return Normal;
1298 }
1299 /*!Returns the barycenter for each element in the support. The barycenter positions are returned
1300   as a field with a number of components equal to the mesh dimension.
1301   The barycenter computed by this method is actually the barycenter of the set of nodes of the elements, each having the same weight.
1302 */
1303 FIELD<double, FullInterlace>* MESH::getBarycenter(const SUPPORT * Support) const throw (MEDEXCEPTION)
1304 {
1305   const char * LOC = "MESH::getBarycenter(SUPPORT*) : ";
1306   const GMESH* myMesh = Support->getMesh();
1307   if (this != myMesh)
1308     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1309
1310   int dim_space = getSpaceDimension();
1311   medEntityMesh support_entity = Support->getEntity();
1312   bool onAll = Support->isOnAllElements();
1313
1314   int nb_type, length_values;
1315   const medGeometryElement* types;
1316   int nb_entity_type;
1317   const int* global_connectivity;
1318
1319   nb_type = Support->getNumberOfTypes();
1320   length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1321   types = Support->getTypes();
1322
1323   FIELD<double, FullInterlace>* Barycenter;
1324   Barycenter = new FIELD<double, FullInterlace>(Support,dim_space);
1325   Barycenter->setName("BARYCENTER");
1326   Barycenter->setDescription("cells or faces barycenter");
1327
1328   for (int k=0;k<dim_space;k++) {
1329     int kp1 = k + 1;
1330     Barycenter->setComponentName(kp1,myMesh->getCoordinatesNames()[k]);
1331     Barycenter->setComponentDescription(kp1,"desc-comp");
1332     Barycenter->setMEDComponentUnit(kp1,myMesh->getCoordinatesUnits()[k]);
1333   }
1334   Barycenter->setIterationNumber(0);
1335   Barycenter->setOrderNumber(0);
1336   Barycenter->setTime(0.0);
1337   double *barycenter=(double *)Barycenter->getValue();
1338   const double * coord = getCoordinates(MED_FULL_INTERLACE);
1339   int index=0;
1340   for (int i=0;i<nb_type;i++)
1341   {
1342     medGeometryElement type = types[i] ;
1343     nb_entity_type = Support->getNumberOfElements(type);
1344     if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA )
1345     {
1346       if (onAll)
1347       {
1348         global_connectivity = getConnectivity(MED_NODAL,support_entity,type);
1349       }
1350       else
1351       {
1352         const int * supp_number = Support->getNumber(type);
1353         const int * connectivity = getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1354         int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1355
1356         for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1357           for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1358             const int *global_connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1359             global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[global_connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1360           }
1361         }
1362         global_connectivity = global_connectivity_tmp;
1363       }
1364     }
1365
1366     switch (type)
1367     {
1368     case MED_TETRA4 : case MED_TETRA10 :
1369       {
1370         for (int tetra =0;tetra<nb_entity_type;tetra++)
1371         {
1372           int tetra_index = (type%100)*tetra;
1373
1374           int N1 = global_connectivity[tetra_index]-1;
1375           int N2 = global_connectivity[tetra_index+1]-1;
1376           int N3 = global_connectivity[tetra_index+2]-1;
1377           int N4 = global_connectivity[tetra_index+3]-1;
1378           double *pts[4];
1379           pts[0]=(double *)coord+dim_space*N1;
1380           pts[1]=(double *)coord+dim_space*N2;
1381           pts[2]=(double *)coord+dim_space*N3;
1382           pts[3]=(double *)coord+dim_space*N4;
1383           INTERP_KERNEL::calculateBarycenter<4,3>((const double **)pts,barycenter+3*index);
1384           index++;
1385         }
1386         break;
1387       }
1388     case MED_PYRA5 : case MED_PYRA13 :
1389       {
1390         for (int pyra=0;pyra<nb_entity_type;pyra++)
1391         {
1392           int pyra_index = (type%100)*pyra;
1393
1394           int N1 = global_connectivity[pyra_index]-1;
1395           int N2 = global_connectivity[pyra_index+1]-1;
1396           int N3 = global_connectivity[pyra_index+2]-1;
1397           int N4 = global_connectivity[pyra_index+3]-1;
1398           int N5 = global_connectivity[pyra_index+4]-1;
1399           double *pts[5];
1400           pts[0]=(double *)coord+dim_space*N1;
1401           pts[1]=(double *)coord+dim_space*N2;
1402           pts[2]=(double *)coord+dim_space*N3;
1403           pts[3]=(double *)coord+dim_space*N4;
1404           pts[4]=(double *)coord+dim_space*N5;
1405           INTERP_KERNEL::calculateBarycenter<5,3>((const double **)pts,barycenter+3*index);
1406           index++;
1407         }
1408         break;
1409       }
1410     case MED_PENTA6 : case MED_PENTA15 :
1411       {
1412         for (int penta=0;penta<nb_entity_type;penta++)
1413         {
1414           int penta_index = (type%100)*penta;
1415
1416           int N1 = global_connectivity[penta_index]-1;
1417           int N2 = global_connectivity[penta_index+1]-1;
1418           int N3 = global_connectivity[penta_index+2]-1;
1419           int N4 = global_connectivity[penta_index+3]-1;
1420           int N5 = global_connectivity[penta_index+4]-1;
1421           int N6 = global_connectivity[penta_index+5]-1;
1422           double *pts[6];
1423           pts[0]=(double *)coord+dim_space*N1;
1424           pts[1]=(double *)coord+dim_space*N2;
1425           pts[2]=(double *)coord+dim_space*N3;
1426           pts[3]=(double *)coord+dim_space*N4;
1427           pts[4]=(double *)coord+dim_space*N5;
1428           pts[5]=(double *)coord+dim_space*N6;
1429           INTERP_KERNEL::calculateBarycenter<6,3>((const double **)pts,barycenter+3*index);
1430           index++;
1431         }
1432         break;
1433       }
1434     case MED_HEXA8 : case MED_HEXA20 :
1435       {
1436         for (int hexa=0;hexa<nb_entity_type;hexa++)
1437         {
1438           int hexa_index = (type%100)*hexa;
1439
1440           int N1 = global_connectivity[hexa_index]-1;
1441           int N2 = global_connectivity[hexa_index+1]-1;
1442           int N3 = global_connectivity[hexa_index+2]-1;
1443           int N4 = global_connectivity[hexa_index+3]-1;
1444           int N5 = global_connectivity[hexa_index+4]-1;
1445           int N6 = global_connectivity[hexa_index+5]-1;
1446           int N7 = global_connectivity[hexa_index+6]-1;
1447           int N8 = global_connectivity[hexa_index+7]-1;
1448           double *pts[8];
1449           pts[0]=(double *)coord+dim_space*N1;
1450           pts[1]=(double *)coord+dim_space*N2;
1451           pts[2]=(double *)coord+dim_space*N3;
1452           pts[3]=(double *)coord+dim_space*N4;
1453           pts[4]=(double *)coord+dim_space*N5;
1454           pts[5]=(double *)coord+dim_space*N6;
1455           pts[6]=(double *)coord+dim_space*N7;
1456           pts[7]=(double *)coord+dim_space*N8;
1457           INTERP_KERNEL::calculateBarycenter<8,3>((const double **)pts,barycenter+3*index);
1458           index++;
1459         }
1460         break;
1461       }
1462     case MED_TRIA3 : case MED_TRIA6 :
1463       {
1464         for (int tria=0;tria<nb_entity_type;tria++)
1465         {
1466           int tria_index = (type%100)*tria;
1467           int N1 = global_connectivity[tria_index]-1;
1468           int N2 = global_connectivity[tria_index+1]-1;
1469           int N3 = global_connectivity[tria_index+2]-1;
1470           double *pts[3];
1471           pts[0]=(double *)coord+dim_space*N1;
1472           pts[1]=(double *)coord+dim_space*N2;
1473           pts[2]=(double *)coord+dim_space*N3;
1474           if (dim_space==2)
1475             INTERP_KERNEL::calculateBarycenter<3,2>((const double **)pts,barycenter+2*index);
1476           else
1477             INTERP_KERNEL::calculateBarycenter<3,3>((const double **)pts,barycenter+3*index);
1478           index++;
1479         }
1480         break;
1481       }
1482     case MED_QUAD4 : case MED_QUAD8 :
1483       {
1484         for (int quad=0;quad<nb_entity_type;quad++)
1485         {
1486           int quad_index = (type%100)*quad;
1487           int N1 = global_connectivity[quad_index]-1;
1488           int N2 = global_connectivity[quad_index+1]-1;
1489           int N3 = global_connectivity[quad_index+2]-1;
1490           int N4 = global_connectivity[quad_index+3]-1;
1491           double *pts[4];
1492           pts[0]=(double *)coord+dim_space*N1;
1493           pts[1]=(double *)coord+dim_space*N2;
1494           pts[2]=(double *)coord+dim_space*N3;
1495           pts[3]=(double *)coord+dim_space*N4;
1496           if (dim_space==2)
1497             INTERP_KERNEL::calculateBarycenter<4,2>((const double **)pts,barycenter+2*index);
1498           else
1499             INTERP_KERNEL::calculateBarycenter<4,3>((const double **)pts,barycenter+3*index);
1500           index++;
1501         }
1502         break;
1503       }
1504     case MED_SEG2 : case MED_SEG3 :
1505       {
1506         for (int edge=0;edge<nb_entity_type;edge++)
1507         {
1508           int edge_index = (type%100)*edge;
1509           int N1 = global_connectivity[edge_index]-1;
1510           int N2 = global_connectivity[edge_index+1]-1;
1511           double *pts[2];
1512           pts[0]=(double *)coord+dim_space*N1;
1513           pts[1]=(double *)coord+dim_space*N2;
1514           if (dim_space==2)
1515             INTERP_KERNEL::calculateBarycenter<2,2>((const double **)pts,barycenter+2*index);
1516           else
1517             INTERP_KERNEL::calculateBarycenter<2,3>((const double **)pts,barycenter+3*index);
1518           index++;
1519         }
1520         break;
1521       }
1522     case MED_POLYGON :
1523       {
1524         if(onAll)
1525         {
1526           int offsetWithClassicType=getGlobalNumberingIndex(support_entity)[i]-1;
1527           const int * connectivity = _connectivity->getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1528           const int * connectivity_index = _connectivity->getConnectivityIndex(MED_EN::MED_NODAL,support_entity) + offsetWithClassicType;
1529           for (int polygs=0;polygs<nb_entity_type;polygs++)
1530           {
1531             int size=connectivity_index[polygs+1]-connectivity_index[polygs];
1532             double **pts=new double * [size];
1533             for(int iPts=0;iPts<size;iPts++)
1534               pts[iPts]=(double *)coord+dim_space*(connectivity[connectivity_index[polygs]+iPts-1]-1);
1535             INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,size,dim_space,barycenter+dim_space*index);
1536             delete [] pts;
1537             index++;
1538           }
1539         }
1540         else
1541         {
1542           const int * supp_number = Support->getNumber(MED_EN::MED_POLYGON);
1543           const int * connectivity = _connectivity->getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1544           const int * connectivity_index = _connectivity->getConnectivityIndex(MED_EN::MED_NODAL,support_entity);
1545           for (int polygs=0;polygs<nb_entity_type;polygs++)
1546           {
1547             int localPolygsNbP1=supp_number[polygs];
1548             int size=connectivity_index[localPolygsNbP1]-connectivity_index[localPolygsNbP1-1];
1549             double **pts=new double * [size];
1550             for(int iPts=0;iPts<size;iPts++)
1551               pts[iPts]=(double *)coord+dim_space*(connectivity[connectivity_index[localPolygsNbP1-1]+iPts-1]-1);
1552             INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,size,dim_space,barycenter+dim_space*index);
1553             delete [] pts;
1554             index++;
1555           }
1556         }
1557         break;
1558       }
1559     case MED_EN::MED_POLYHEDRA:
1560       {
1561         if(onAll)
1562         {
1563           for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
1564           {
1565             int lgthNodes;
1566             int offsetWithClassicType=getGlobalNumberingIndex(support_entity)[i]-1;
1567             int *nodes=_connectivity->getNodesOfPolyhedron(offsetWithClassicType+polyhs+1,lgthNodes);
1568             double **pts=new double * [lgthNodes];
1569             for(int iPts=0;iPts<lgthNodes;iPts++)
1570               pts[iPts]=(double *)coord+3*(nodes[iPts]-1);
1571             INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,lgthNodes,3,barycenter+3*index);
1572             delete [] pts;
1573             delete [] nodes;
1574             index++;
1575           }
1576         }
1577         else
1578         {
1579           const int * supp_number = Support->getNumber(MED_EN::MED_POLYHEDRA);
1580           for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
1581           {
1582             int lgthNodes;
1583             int *nodes=_connectivity->getNodesOfPolyhedron(supp_number[polyhs],lgthNodes);
1584             double **pts=new double * [lgthNodes];
1585             for(int iPts=0;iPts<lgthNodes;iPts++)
1586               pts[iPts]=(double *)coord+3*(nodes[iPts]-1);
1587             INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,lgthNodes,3,barycenter+3*index);
1588             delete [] pts;
1589             delete [] nodes;
1590             index++;
1591           }
1592         }
1593         break;
1594       }
1595     default :
1596       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get a barycenter on it (in fact unknown type) !"));
1597       break;
1598     }
1599
1600     if (!onAll)
1601       if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
1602         delete [] global_connectivity;
1603   }
1604   return Barycenter;
1605 }
1606 /*!
1607   @}
1608 */
1609
1610 bool MESH::isEmpty() const
1611 {
1612   bool notempty = _name != "NOT DEFINED"                || _coordinate != NULL           || _connectivity != NULL ||
1613     _spaceDimension !=MED_INVALID ||
1614     _numberOfNodes !=MED_INVALID  || _groupNode.size() != 0   ||
1615     _familyNode.size() != 0       || _groupCell.size() != 0   ||
1616     _familyCell.size() != 0       || _groupFace.size() != 0   ||
1617     _familyFace.size() != 0       || _groupEdge.size() != 0   ||
1618     _familyEdge.size() != 0;
1619   return !notempty;
1620 }
1621
1622 //================================================================================
1623 /*!
1624  * \brief Check nature of GMESH
1625  */
1626 //================================================================================
1627
1628 bool MESH::getIsAGrid() const
1629 {
1630   return false;
1631 }
1632
1633 /*!
1634  * \brief Implement pure virtual method used to get MESH from GMESH
1635  */
1636 const MESH* MESH::convertInMESH() const
1637 {
1638   this->addReference();
1639   return this;
1640 }
1641
1642 /*!
1643 \addtogroup MESH_advanced
1644 @{
1645 */
1646
1647 /*!
1648   Retrieves the skin of support \a Support3D. This method is only available in 3D.
1649   On output, it returns a MED_FACE support with the skin of all elements contained in support.
1650   The skin is defined as the list of faces that are compnents of only one volume in the input
1651   support.
1652
1653   WARNING: This method can recalculate descending connectivity from partial to full form,
1654   so that partial SUPPORT on MED_FACE on this mesh becomes invalid.
1655 */
1656 SUPPORT * MESH::getSkin(const SUPPORT * Support3D) throw (MEDEXCEPTION)
1657 {
1658   const char * LOC = "MESH::getSkin : " ;
1659   BEGIN_OF_MED(LOC) ;
1660   // some test :
1661   if (this != Support3D->getMesh())
1662     throw MEDEXCEPTION(STRING(LOC) <<  "no compatibility between *this and SUPPORT::_mesh !");
1663   if (getMeshDimension() != 3 || Support3D->getEntity() != MED_CELL)
1664     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Defined on 3D cells only"));
1665
1666   // well, all rigth !
1667   SUPPORT * mySupport = new SUPPORT;
1668   mySupport->setMesh((MESH *)this);
1669   mySupport->setName("Skin");
1670   mySupport->setEntity( MED_FACE );
1671
1672   list<int> myElementsList;
1673   int i,j, size = 0;
1674
1675   // assure that descending connectivity is full
1676   if ( !_connectivity )
1677     throw MEDEXCEPTION(STRING(LOC) << "no connectivity defined in MESH !");
1678   _connectivity->calculateFullDescendingConnectivity(MED_CELL);
1679   if (Support3D->isOnAllElements())
1680   {
1681     const int* value = getReverseConnectivity(MED_DESCENDING);
1682     const int* index = getReverseConnectivityIndex(MED_DESCENDING);
1683
1684     int nbFaces = getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS);
1685     for (int i = 0; i < nbFaces; i++)
1686     {
1687       //a face is in skin if it is linked to one element or if one of the elements
1688       //it is linked to is "virtual"
1689       if ((index[i+1]-index[i]==1) || (value[index[i]-1]==0) || (value[index[i]]==0)) {
1690         myElementsList.push_back( i+1 );
1691         size++;
1692       }
1693     }
1694   }
1695   else
1696   {
1697     map<int,int> FaceNbEncounterNb;
1698
1699     int * myConnectivityValue = const_cast <int *>
1700       (getConnectivity(MED_DESCENDING,MED_CELL, MED_ALL_ELEMENTS));
1701     int * myConnectivityIndex = const_cast <int *> (getConnectivityIndex(MED_DESCENDING, MED_CELL));
1702     int * myCellNbs = const_cast <int *> (Support3D->getnumber()->getValue());
1703     int nbCells = Support3D->getnumber()->getLength();
1704     for (i=0; i<nbCells; ++i)
1705     {
1706       int cellNb = myCellNbs [ i ];
1707       int faceFirst = myConnectivityIndex[ cellNb-1 ];
1708       int faceLast  = myConnectivityIndex[ cellNb ];
1709       for (j = faceFirst; j < faceLast; ++j)
1710       {
1711         int faceNb = abs( myConnectivityValue [ j-1 ] );
1712         if (FaceNbEncounterNb.find( faceNb ) == FaceNbEncounterNb.end())
1713           FaceNbEncounterNb[ faceNb ] = 1;
1714         else
1715           FaceNbEncounterNb[ faceNb ] ++;
1716       }
1717     }
1718     map<int,int>::iterator FaceNbEncounterNbItr;
1719     for (FaceNbEncounterNbItr = FaceNbEncounterNb.begin();
1720          FaceNbEncounterNbItr != FaceNbEncounterNb.end();
1721          FaceNbEncounterNbItr ++)
1722       if ((*FaceNbEncounterNbItr).second == 1)
1723       {
1724         myElementsList.push_back( (*FaceNbEncounterNbItr).first) ;
1725         size++ ;
1726       }
1727   }
1728   // Well, we must know how many geometric type we have found
1729   int * myListArray = new int[size] ;
1730   int id = 0 ;
1731   list<int>::iterator myElementsListIt ;
1732   for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
1733     myListArray[id]=(*myElementsListIt) ;
1734     id ++ ;
1735   }
1736
1737   int numberOfGeometricType ;
1738   medGeometryElement* geometricType ;
1739   int * geometricTypeNumber ;
1740   int * numberOfEntities ;
1741   int * mySkyLineArrayIndex ;
1742
1743   int numberOfType = getNumberOfTypes(MED_FACE) ;
1744   if (numberOfType == 1) // wonderfull : it's easy !
1745     {
1746       numberOfGeometricType = 1 ;
1747       geometricType = new medGeometryElement[1] ;
1748       const medGeometryElement *  allType = getTypes(MED_FACE);
1749       geometricType[0] = allType[0] ;
1750       geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
1751       geometricTypeNumber[0] = 0 ;
1752       numberOfEntities = new int[1] ;
1753       numberOfEntities[0] = size ;
1754       mySkyLineArrayIndex = new int[2] ;
1755       mySkyLineArrayIndex[0]=1 ;
1756       mySkyLineArrayIndex[1]=1+size ;
1757     }
1758   else // hemmm
1759     {
1760       map<medGeometryElement,int> theType ;
1761       for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++)
1762         {
1763           medGeometryElement myType = getElementType(MED_FACE,*myElementsListIt) ;
1764           if (theType.find(myType) != theType.end() )
1765             theType[myType]+=1 ;
1766           else
1767             theType[myType]=1 ;
1768         }
1769       numberOfGeometricType = theType.size() ;
1770       geometricType = new medGeometryElement[numberOfGeometricType] ;
1771       geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
1772       numberOfEntities = new int[numberOfGeometricType] ;
1773       mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
1774       int index = 0 ;
1775       mySkyLineArrayIndex[0]=1 ;
1776       map<medGeometryElement,int>::iterator theTypeIt ;
1777       for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++)
1778         {
1779           geometricType[index] = (*theTypeIt).first ;
1780           geometricTypeNumber[index] = 0 ;
1781           numberOfEntities[index] = (*theTypeIt).second ;
1782           mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfEntities[index] ;
1783           index++ ;
1784         }
1785     }
1786   MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
1787
1788   mySupport->setNumberOfGeometricType(numberOfGeometricType) ;
1789   mySupport->setGeometricType(geometricType) ;
1790   mySupport->setNumberOfElements(numberOfEntities) ;
1791   mySupport->setNumber(mySkyLineArray) ;
1792
1793   delete[] numberOfEntities;
1794   delete[] geometricTypeNumber;
1795   delete[] geometricType;
1796   delete[] mySkyLineArrayIndex;
1797   delete[] myListArray;
1798
1799   END_OF_MED(LOC);
1800   return mySupport ;
1801
1802 }
1803
1804
1805
1806 int MESH::getElementContainingPoint(const double *coord)
1807 {
1808   if(_spaceDimension!=3 && _spaceDimension!=2)
1809     throw MEDEXCEPTION("MESH::getElementContainingPoint : invalid _spaceDimension must be equal to 2 or 3 !!!");
1810   PointLocator loc(*this);
1811   std::list<int> li=loc.locate(coord);
1812   if(li.empty())
1813     return -1;
1814   return li.front();
1815 }
1816
1817 //! Converts MED_CELL connectivity to polyhedra connectivity
1818 //! Converts MED_FACE connectivity to polygon connectivity
1819 //! Wil work only for 3D meshes with nodal connectivity
1820
1821 void MESH::convertToPoly()
1822 {
1823   if (getMeshDimension()!=3) return;
1824
1825   CONNECTIVITY* newpolygonconnectivity = 0;
1826   CONNECTIVITY* newpolyhedraconnectivity = new CONNECTIVITY(MED_EN::MED_CELL);
1827
1828   {
1829     ////////////////////////////////////////////:
1830     // First step : Treating polygons connectivity
1831     ///////////////////////////////////////////
1832
1833     int oldnbface = getNumberOfElements(MED_EN::MED_FACE,MED_EN::MED_ALL_ELEMENTS);
1834     int nbTypes = oldnbface > 0 ? 1 : 0;
1835     newpolygonconnectivity = new CONNECTIVITY(nbTypes, MED_EN::MED_FACE);
1836     if ( nbTypes > 0 )
1837     {
1838       medGeometryElement type = MED_POLYGON;
1839       newpolygonconnectivity->setGeometricTypes( &type, MED_FACE );
1840
1841       const int count[] = { 1, oldnbface + 1 };
1842       newpolygonconnectivity->setCount( count, MED_FACE );
1843
1844       const int* oldconn = getConnectivity(MED_EN::MED_NODAL, MED_EN::MED_FACE, MED_EN::MED_ALL_ELEMENTS);
1845       const int* oldconnindex= getConnectivityIndex(MED_EN::MED_NODAL,MED_EN::MED_FACE);
1846       newpolygonconnectivity->setNodal( oldconn, MED_FACE, type, oldconnindex );
1847
1848       newpolygonconnectivity->setNumberOfNodes( getNumberOfNodes() );
1849       newpolygonconnectivity->setEntityDimension( 2 );
1850     }
1851   }
1852   ///////////////////////////////////////////
1853   // 2nd step : Treating polyhedra connectivity
1854   //////////////////////////////////////////
1855   {
1856     vector<int> newconn;
1857     vector<int> newindex(1,1);
1858
1859     int nboldtypes=getNumberOfTypes(MED_EN::MED_CELL);
1860     const MED_EN::medGeometryElement* oldtypes = getTypes(MED_EN::MED_CELL);
1861     const int* oldconn = getConnectivity(MED_EN::MED_NODAL, MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS);
1862     const int* oldconnindex= getConnectivityIndex(MED_EN::MED_NODAL,MED_EN::MED_CELL);
1863     int oldnbelem = getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
1864
1865     for (int itype=0; itype<nboldtypes; itype++)
1866     {
1867       MED_EN::medGeometryElement type = oldtypes[itype];
1868       int nb_elems=getNumberOfElements(MED_EN::MED_CELL,type);
1869       if ( type == MED_POLYHEDRA )
1870       {
1871         const int* oldpolyindex = oldconnindex + getGlobalNumberingIndex( MED_CELL )[itype] - 1;
1872         int oldpolyconnsize = oldpolyindex[nb_elems] - oldpolyindex[0];
1873         newconn.insert( newconn.end(), oldconn, oldconn + oldpolyconnsize );
1874         int delta = newindex.back() - oldpolyindex[0];
1875         for (int ielem=0; ielem<nb_elems; ielem++)
1876           newindex.push_back( delta + oldpolyindex[ ielem+1 ]);
1877       }
1878       else
1879       {
1880         MEDMEM::CELLMODEL cellmodel(type);
1881         int nbfacespertype = cellmodel.getNumberOfConstituents(1);
1882         int nbnodespertype = cellmodel.getNumberOfNodes();
1883         for (int ielem=0; ielem<nb_elems; ielem++)
1884         {
1885           for (int iface=0; iface< nbfacespertype; iface++)
1886           {
1887             //local conn contains the local nodal connectivity for the iface-th face of type
1888             const int* local_conn = cellmodel.getNodesConstituent(1,iface+1);
1889             medGeometryElement facetype = cellmodel.getConstituentType(1,iface+1);
1890             int nbface_nodes=facetype%100;
1891             for ( int inode=0; inode<nbface_nodes;inode++)
1892               newconn.push_back(oldconn[local_conn[inode]-1]);
1893             if ( iface != nbfacespertype-1 )
1894               newconn.push_back(-1);
1895           }
1896           newindex.push_back( newconn.size() + 1 );
1897           oldconn += nbnodespertype;
1898         }
1899       }
1900     }
1901     int nbTypes = oldnbelem > 0 ? 1 : 0;
1902     if ( newpolyhedraconnectivity ) delete newpolyhedraconnectivity;
1903     newpolyhedraconnectivity = new CONNECTIVITY(nbTypes, MED_EN::MED_CELL);
1904     if ( nbTypes > 0 )
1905     {
1906       medGeometryElement type = MED_POLYHEDRA;
1907       newpolyhedraconnectivity->setGeometricTypes( &type, MED_CELL );
1908
1909       const int count[] = { 1, oldnbelem + 1 };
1910       newpolyhedraconnectivity->setCount( count, MED_CELL );
1911
1912       newpolyhedraconnectivity->setNodal( &newconn[0], MED_CELL, type, &newindex[0] );
1913
1914       newpolyhedraconnectivity->setNumberOfNodes( getNumberOfNodes() );
1915       newpolyhedraconnectivity->setEntityDimension( 3 );
1916     }
1917   }
1918
1919   delete _connectivity;
1920
1921   _connectivity=newpolyhedraconnectivity;
1922   _connectivity->setConstituent(newpolygonconnectivity);
1923
1924 }
1925
1926 vector< vector<double> > MESH::getBoundingBox() const
1927 {
1928   const double *myCoords=_coordinate->getCoordinates(MED_EN::MED_FULL_INTERLACE);
1929   vector< vector<double> > ret(2);
1930   int i,j;
1931   ret[0].resize(_spaceDimension);
1932   ret[1].resize(_spaceDimension);
1933   for(i=0;i<_spaceDimension;i++)
1934   {
1935     ret[0][i]=1.e300;
1936     ret[1][i]=-1.e300;
1937   }
1938   for(i=0;i<_coordinate->getNumberOfNodes();i++)
1939   {
1940     for(j=0;j<_spaceDimension;j++)
1941     {
1942       double tmp=myCoords[i*_spaceDimension+j];
1943       if(tmp<ret[0][j])
1944         ret[0][j]=tmp;
1945       if(tmp>ret[1][j])
1946         ret[1][j]=tmp;
1947     }
1948   }
1949   return ret;
1950 }