src/MEDMEM/MEDMEM_Grid.cxx

   1 // Copyright (C) 2007-2013  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 // File      : MEDMEM_Grid.hxx
  24 // Created   : Wed Dec 18 08:35:26 2002
  25 // Descr     : class containing structured mesh data
  26 // Author    : Edward AGAPOV (eap)
  27 // Project   : SALOME Pro
  28 // Module    : MED
  29 //
  30 #include "MEDMEM_Grid.hxx"
  31 #include "MEDMEM_Meshing.hxx"
  32 #include "MEDMEM_CellModel.hxx"
  33 #include "MEDMEM_SkyLineArray.hxx"
  34 #include "MEDMEM_DriverFactory.hxx"
  35
  36 #include <memory>
  37
  38 using namespace std;
  39 using namespace MEDMEM;
  40 using namespace MED_EN;
  41
  42
  43 // Block defining the comments for the MEDMEM_ug documentation
  44
  45 /*!
  46
  47 \defgroup GRID_axes Information about axes
  48 This group of methods retrieves information about the axes of the grid.
  49
  50 \defgroup GRID_connectivity Utility methods for defining element positions in the grid
  51 These methods enable the user to convert a position on the grid to a global element number
  52
  53 \defgroup GRID_constructors Constructors
  54 These methods are the different constructors for the grid objects.
  55
  56 */
  57
  58 namespace
  59 {
  60   const string* defaultStrings()
  61   {
  62     static const string str[3] = { "UNDEFINED", "UNDEFINED", "UNDEFINED" };
  63     return str;
  64   }
  65 }
  66
  67 //=======================================================================
  68 //function : GRID
  69 //purpose  : empty constructor
  70 //=======================================================================
  71
  72 GRID::GRID() {
  73   init();
  74   MESSAGE_MED("A GRID CREATED");
  75 }
  76 //
  77 //=======================================================================
  78 ////function : GRID
  79 ////purpose  : array constructor
  80 ////=======================================================================
  81 /*!
  82 \if MEDMEM_ug
  83 \addtogroup GRID_constructors
  84 @{
  85 \endif
  86 */
  87 /*!
  88  * \brief Constructor specifying the axes of the grid
  89  *
  90  * This constructor describes the grid by specifying the location of the nodes on
  91 each of the axis. The dimension of the grid is implicitly defined by the
  92 size of vector \a xyz_array.
  93  *
  94  *\param xyz_array specifies the node coordinates for each direction
  95  *\param coord_name names of the different coordinates
  96  *\param coord_unit names of the different coordinate units
  97  *\param type  grid type (MED_POLAR, MED_CARTESIAN)
  98 */
  99 GRID::GRID(const std::vector<std::vector<double> >& xyz_array,
 100            const std::vector<std::string>&          coord_name,
 101            const std::vector<std::string>&          coord_unit,
 102            const MED_EN::med_grid_type              type)
 103   :_gridType(type)
 104 {
 105     init(); // PAL 12136
 106     _is_default_gridType = false;
 107
 108     _spaceDimension = xyz_array.size();
 109
 110     // create a non allocated COORDINATE
 111     _coordinate = new COORDINATE(_spaceDimension, &coord_name[0], &coord_unit[0]);
 112     string coordinateSystem = "UNDEFINED";
 113     if( _gridType == MED_CARTESIAN)
 114       coordinateSystem = "CARTESIAN";
 115     else if ( _gridType == MED_POLAR)
 116       coordinateSystem = "CYLINDRICAL";
 117     _coordinate->setCoordinatesSystem(coordinateSystem);
 118
 119     // set the GRID part
 120     if (_spaceDimension>=1)
 121     {
 122         _iArrayLength=xyz_array[0].size();
 123         _iArray=new double[_iArrayLength];
 124         std::copy(xyz_array[0].begin(),xyz_array[0].end(),_iArray);
 125     }
 126     if (_spaceDimension>=2)
 127     {
 128         _jArrayLength=xyz_array[1].size();
 129         _jArray=new double[_jArrayLength];
 130         std::copy(xyz_array[1].begin(),xyz_array[1].end(),_jArray);
 131     }
 132     if (_spaceDimension>=3)
 133     {
 134         _kArrayLength=xyz_array[2].size();
 135         _kArray=new double[_kArrayLength];
 136         std::copy(xyz_array[2].begin(),xyz_array[2].end(),_kArray);
 137     }
 138 }
 139
 140 //================================================================================
 141 /*!
 142  * \brief Reads a GRID form the file
 143  *  \param driverType - type of driver to read the specified file
 144  *  \param fileName - the file name to read
 145  *  \param driverName - name of a grid to read
 146  */
 147 //================================================================================
 148
 149 GRID::GRID(driverTypes driverType, const string &  fileName, const string &  driverName)
 150 {
 151
 152   const char* LOC = "GRID::GRID(driverTypes , const string &  , const string &) : ";
 153   BEGIN_OF_MED(LOC);
 154
 155   init();
 156   auto_ptr<GENDRIVER> myDriver (DRIVERFACTORY::buildDriverForMesh(driverType,fileName,this,driverName,RDONLY));
 157   myDriver->open();
 158   myDriver->read();
 159   myDriver->close();
 160
 161   END_OF_MED(LOC);
 162 }
 163
 164 /*!\if MEDMEM_ug @} \endif */
 165
 166 //=======================================================================
 167 //function : GRID
 168 //purpose  : empty constructor
 169 //=======================================================================
 170
 171 GRID::GRID(const MED_EN::med_grid_type type)
 172 {
 173   init();
 174   _gridType = type;
 175   _is_default_gridType = false;
 176   MESSAGE_MED("A TYPED GRID CREATED");
 177 }
 178
 179 //=======================================================================
 180 //function : GRID
 181 //purpose  : copying constructor
 182 //=======================================================================
 183
 184 GRID::GRID(const GRID& otherGrid) {
 185   *this = otherGrid;
 186 }
 187
 188 //=======================================================================
 189 //function : ~GRID
 190 //purpose  : destructor
 191 //=======================================================================
 192
 193 GRID::~GRID() {
 194   MESSAGE_MED("GRID::~GRID() : Destroying the Grid");
 195   if ( _coordinate ) delete _coordinate; _coordinate = 0;
 196   if ( _iArray != (double* ) NULL) delete [] _iArray;
 197   if ( _jArray != (double* ) NULL) delete [] _jArray;
 198   if ( _kArray != (double* ) NULL) delete [] _kArray;
 199 }
 200
 201 //=======================================================================
 202 //function : init
 203 //purpose :
 204 //=======================================================================
 205
 206 void GRID::init()
 207 {
 208   GMESH::init();
 209
 210   _gridType = MED_CARTESIAN;
 211   _is_default_gridType = true;
 212   _coordinate = 0;
 213   _iArray = _jArray = _kArray = (double* ) NULL;
 214   _iArrayLength = _jArrayLength = _kArrayLength = 0;
 215
 216 }
 217
 218 //================================================================================
 219 /*!
 220  * \brief Return true if contains no elements
 221  */
 222 //================================================================================
 223
 224 bool GRID::isEmpty() const
 225 {
 226   return !_iArrayLength && !_coordinate;
 227 }
 228
 229
 230 //=======================================================================
 231 //function: operator=
 232 //purpose : operator=
 233 //=======================================================================
 234
 235 GRID & GRID::operator=(const GRID & otherGrid)
 236 {
 237   // NOT IMPLEMENTED
 238
 239   GMESH::operator=(otherGrid);
 240   return *this;
 241 }
 242
 243 //=======================================================================
 244 /*!
 245   Returns true if mesh \a other has same
 246   coordinates (to 1E-15 precision ) and same connectivity as the calling object.
 247   Information like name or description is not taken into account
 248   for the comparison.
 249 */
 250 //=======================================================================
 251
 252 bool GRID::deepCompare(const GMESH& gother) const
 253 {
 254   if ( gother.getIsAGrid() != getIsAGrid())
 255     return false;
 256   const GRID& other = static_cast<const GRID&>( gother );
 257
 258   if ( getSpaceDimension() != other.getSpaceDimension() )
 259     return false;
 260
 261   if ( _gridType != other._gridType )
 262     return false;
 263
 264   if( bool( _coordinate) != bool(other._coordinate))
 265     return false;
 266
 267   if ( _coordinate->getNumberOfNodes() > 0 )
 268   {
 269     if ( _coordinate->getNumberOfNodes() != other._coordinate->getNumberOfNodes() )
 270       return false;
 271     const int size = _coordinate->getNumberOfNodes() * getSpaceDimension();
 272     const double* coord1=_coordinate->getCoordinates(MED_FULL_INTERLACE);
 273     const double* coord2=other._coordinate->getCoordinates(MED_FULL_INTERLACE);
 274     for(int i = 0 ; i < size; i++ )
 275       if ( !(fabs(coord1[i]-coord2[i])<1e-15))
 276         return false;
 277   }
 278
 279   if ( _iArrayLength != other._iArrayLength )
 280     return false;
 281   if ( _jArrayLength != other._jArrayLength )
 282     return false;
 283   if ( _kArrayLength != other._kArrayLength )
 284     return false;
 285
 286   if ( bool(_iArray) != bool(other._iArray) )
 287     return false;
 288   if ( bool(_jArray) != bool(other._jArray) )
 289     return false;
 290   if ( bool(_kArray) != bool(other._kArray) )
 291     return false;
 292
 293   if ( _iArray )
 294     for ( int i = 0; i < _iArrayLength; ++i )
 295       if ( !(fabs(_iArray[i]-other._iArray[i])<1e-15))
 296         return false;
 297   if ( _jArray )
 298     for ( int i = 0; i < _jArrayLength; ++i )
 299       if ( !(fabs(_jArray[i]-other._jArray[i])<1e-15))
 300         return false;
 301   if ( _kArray )
 302     for ( int i = 0; i < _kArrayLength; ++i )
 303       if ( !(fabs(_kArray[i]-other._kArray[i])<1e-15))
 304         return false;
 305
 306   return true;
 307 }
 308
 309 //================================================================================
 310 /*!
 311  * \brief print my contents
 312  */
 313 //================================================================================
 314
 315 void GRID::printMySelf(std::ostream &os) const
 316 {
 317   // TODO
 318   cout << "NOT IMPLEMENTED" << endl;
 319 }
 320
 321 //================================================================================
 322 /*!
 323  * \brief Create an unstructured MESH. Call removeReference() after having finished using it!!!
 324  */
 325 //================================================================================
 326
 327 const MESH * GRID::convertInMESH() const
 328 {
 329   MESHING* mesh = new MESHING;
 330   mesh->setName( getName() );
 331
 332   int i, j, k;
 333
 334   // ---------------
 335   // 1. Coordinates
 336   // ---------------
 337
 338   PointerOf< double > coords;
 339   if ( _gridType == MED_BODY_FITTED )
 340   {
 341     if ( !_coordinate )
 342       throw MEDEXCEPTION
 343         (LOCALIZED("GRID::convertInMESH() : _coordinate of MED_BODY_FITTED GRID not defined !"));
 344     coords.set( _coordinate->getCoordinates(MED_FULL_INTERLACE));
 345   }
 346   else
 347   {
 348     coords.set( getSpaceDimension() * getNumberOfNodes() );
 349     double* coord = coords;
 350     switch ( getSpaceDimension() )
 351     {
 352     case 3:
 353       for (k=0; k < _kArrayLength; ++k)
 354         for (j=0; j < _jArrayLength; ++j)
 355           for (i=0; i < _iArrayLength; ++i)
 356             *coord++ = _iArray[i], *coord++ = _jArray[j], *coord++ = _kArray[k];
 357       break;
 358
 359     case 2:
 360       for (j=0; j < _jArrayLength; ++j)
 361         for (i=0; i < _iArrayLength; ++i)
 362           *coord++ = _iArray[i], *coord++ = _jArray[j];
 363       break;
 364
 365     case 1:
 366       coords.set(_iArray);
 367       break;
 368     }
 369   }
 370   mesh->setCoordinates( getSpaceDimension(),
 371                         getNumberOfNodes(),
 372                         (const double *) coords,
 373                         getCoordinatesSystem(),
 374                         MED_EN::MED_FULL_INTERLACE);
 375   mesh->setCoordinatesNames( getCoordinatesNames() );
 376   mesh->setCoordinatesUnits( getCoordinatesUnits() );
 377
 378   // ----------------
 379   // 2. Connectivity
 380   // ----------------
 381
 382   // 2.1 Cells
 383   // ---------
 384   medEntityMesh subEntity;
 385   {
 386     mesh->setNumberOfTypes( getNumberOfTypes(MED_CELL), MED_CELL );
 387     mesh->setTypes( getTypes( MED_CELL), MED_CELL );
 388     int nbCells = getNumberOfElements( MED_CELL, MED_ALL_ELEMENTS );
 389     mesh->setNumberOfElements( &nbCells, MED_CELL );
 390
 391     vector<int> conn;
 392     switch ( _spaceDimension )
 393     {
 394     case 3: // HEXA8
 395       for ( k = 0; k < _kArrayLength-1; ++k )
 396         for ( j = 0; j < _jArrayLength-1; ++j )
 397           for ( i = 0; i < _iArrayLength-1; ++i )
 398           {
 399             conn.push_back( getNodeNumber( i  , j  , k  ));
 400             conn.push_back( getNodeNumber( i  , j+1, k  ));
 401             conn.push_back( getNodeNumber( i+1, j+1, k  ));
 402             conn.push_back( getNodeNumber( i+1, j  , k  ));
 403             conn.push_back( getNodeNumber( i  , j  , k+1));
 404             conn.push_back( getNodeNumber( i  , j+1, k+1));
 405             conn.push_back( getNodeNumber( i+1, j+1, k+1));
 406             conn.push_back( getNodeNumber( i+1, j  , k+1));
 407           }
 408       subEntity = MED_FACE;
 409       break;
 410
 411     case 2: // QUAD4
 412       for ( j = 0; j < _jArrayLength-1; ++j )
 413         for ( i = 0; i < _iArrayLength-1; ++i )
 414         {
 415           int n1 = 1 + i + j*_iArrayLength;
 416           conn.push_back( n1 );
 417           conn.push_back( n1 + _iArrayLength );
 418           conn.push_back( n1 + _iArrayLength + 1 );
 419           conn.push_back( n1 + 1 );
 420         }
 421       subEntity = MED_EDGE;
 422       break;
 423
 424     case 1: // SEG2
 425       for ( i = 0; i < _iArrayLength-1; ++i )
 426       {
 427         conn.push_back( i + 1 );
 428         conn.push_back( i + 2 );
 429       }
 430       break;
 431     }
 432     mesh->setConnectivity( MED_CELL, getTypes(MED_CELL)[0], &conn[0] );
 433   }
 434
 435   // 2.2 sub entities
 436   // -----------------
 437
 438   if ( _spaceDimension > 1 )
 439   {
 440     mesh->setNumberOfTypes( getNumberOfTypes(subEntity), subEntity );
 441     mesh->setTypes( getTypes( subEntity), subEntity );
 442     int nbCells = getNumberOfElements( subEntity, MED_ALL_ELEMENTS );
 443     mesh->setNumberOfElements( &nbCells, subEntity );
 444
 445     vector<int> conn;
 446     switch ( _spaceDimension )
 447     {
 448     case 3: // QUAD4
 449
 450       // normal to OX
 451       for ( k = 0; k < _kArrayLength-1; ++k )
 452         for ( j = 0; j < _jArrayLength-1; ++j )
 453           for ( i = 0; i < _iArrayLength; ++i )
 454           {
 455             conn.push_back( getNodeNumber( i, j  , k  ));
 456             conn.push_back( getNodeNumber( i, j  , k+1));
 457             conn.push_back( getNodeNumber( i, j+1, k+1));
 458             conn.push_back( getNodeNumber( i, j+1, k  ));
 459           }
 460       // normal to OY
 461       for ( k = 0; k < _kArrayLength-1; ++k )
 462         for ( j = 0; j < _jArrayLength; ++j )
 463           for ( i = 0; i < _iArrayLength-1; ++i )
 464           {
 465             conn.push_back( getNodeNumber( i  , j, k  ));
 466             conn.push_back( getNodeNumber( i+1, j, k  ));
 467             conn.push_back( getNodeNumber( i+1, j, k+1));
 468             conn.push_back( getNodeNumber( i  , j, k+1));
 469           }
 470       // normal to OZ
 471       for ( k = 0; k < _kArrayLength; ++k )
 472         for ( j = 0; j < _jArrayLength-1; ++j )
 473           for ( i = 0; i < _iArrayLength-1; ++i )
 474           {
 475             conn.push_back( getNodeNumber( i  , j  , k));
 476             conn.push_back( getNodeNumber( i  , j+1, k));
 477             conn.push_back( getNodeNumber( i+1, j+1, k));
 478             conn.push_back( getNodeNumber( i+1, j  , k));
 479           }
 480       break;
 481
 482     case 2: // SEG2
 483
 484       // || OX
 485       for ( j = 0; j < _jArrayLength; ++j )
 486         for ( i = 0; i < _iArrayLength-1; ++i )
 487         {
 488           int n1 = 1 + i + j*_iArrayLength;
 489           conn.push_back( n1 );
 490           conn.push_back( n1 + 1 );
 491         }
 492       // || OY
 493       for ( j = 0; j < _jArrayLength-1; ++j )
 494         for ( i = 0; i < _iArrayLength; ++i )
 495         {
 496           int n1 = 1 + i + j*_iArrayLength;
 497           conn.push_back( n1 );
 498           conn.push_back( n1 + _iArrayLength );
 499         }
 500       break;
 501
 502     }
 503     mesh->setConnectivity( subEntity, getTypes(subEntity)[0], &conn[0] );
 504   }
 505
 506   // 3. Groups and Families
 507   // -----------------------
 508
 509   const vector<GROUP*>* groups[] = { &_groupNode, &_groupCell, &_groupFace, &_groupEdge };
 510   for ( i = 0; i < 4; ++i )
 511     for ( j = 0; j < (int)groups[i]->size(); ++j )
 512       mesh->addGroup( * groups[i]->at( j ));
 513
 514   return mesh;
 515 }
 516
 517 //=======================================================================
 518 /*!
 519   return the GRID Geometric type, without computing all connectivity
 520 */
 521 //=======================================================================
 522
 523 const medGeometryElement * GRID::getTypes(MED_EN::medEntityMesh entity) const
 524 {
 525   static const medGeometryElement _gridGeometry[4]={MED_HEXA8,MED_QUAD4,MED_SEG2,MED_POINT1};
 526   int i=0;
 527   if(entity==MED_CELL)
 528   {
 529     i=3-_spaceDimension;
 530   }
 531   else if(entity==MED_FACE && _spaceDimension>2 )
 532     i=1;
 533   else if(entity==MED_EDGE && _spaceDimension>1 )
 534     i=2;
 535   else if(entity==MED_NODE && _spaceDimension>0)
 536     i=3;
 537   else
 538     throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::getGeometricTypes : Entity not defined !"));
 539   return &_gridGeometry[i];
 540 }
 541
 542 //=======================================================================
 543 //function : getArrayLength
 544 //purpose  : return array length. Axis = [1,2,3] meaning [i,j,k],
 545 //=======================================================================
 546 /*!\if MEDMEM_ug
 547 \addtogroup GRID_axes
 548 @{
 549 \endif
 550 */
 551 /*! Returns the number of nodes on axis number \a Axis (axis numbering starts at 1).
 552 */
 553 int GRID::getArrayLength( const int Axis ) const throw (MEDEXCEPTION)
 554 {
 555   switch (Axis) {
 556   case 1: return _iArrayLength;
 557   case 2: return _jArrayLength;
 558   case 3: return _kArrayLength;
 559   default:
 560     throw MED_EXCEPTION ( LOCALIZED( STRING("GRID::getArrayLength ( ") << Axis << ")"));
 561   }
 562   return 0;
 563 }
 564
 565 //=======================================================================
 566 //function : getArrayValue
 567 //purpose  : return i-th array component. Axis = [1,2,3] meaning [i,j,k],
 568 //           exception if Axis out of [1-3] range
 569 //           exception if i is out of range 0 <= i < getArrayLength(Axis);
 570 //=======================================================================
 571 /*!
 572 Returns the value of node coordinate \a i on axis \a Axis.
 573  */
 574 const double GRID::getArrayValue (const int Axis, const int i) const throw (MEDEXCEPTION)
 575 {
 576   if (i < 0 || i >= getArrayLength(Axis))
 577     throw MED_EXCEPTION
 578       ( LOCALIZED(STRING("GRID::getArrayValue ( ") << Axis << ", " << i << ")"));
 579
 580   switch (Axis) {
 581   case 1:  return _iArray[ i ];
 582   case 2:  return _jArray[ i ];
 583   default: return _kArray[ i ];
 584   }
 585   return 0.0;
 586 }
 587
 588 /*!
 589 \if MEDMEM_ug
 590 @}
 591 \endif
 592  */
 593
 594 /*!
 595 \if MEDMEM_ug
 596 \addtogroup GRID_connectivity
 597 @{
 598 \endif
 599 */
 600 /*!
 601 @{
 602 \name Position to number conversion methods
 603 \a getXXXNumber methods enable the user to convert an \f$ (i,j,k)\f$ position into a global number in the array.
 604
 605 Axis [1,2,3] means one of directions: along \a i, \a j or \a k.
 606 For cell constituents (FACE or EDGE), Axis selects one of those having same  \f$ (i, j, k )\f$ :
 607 - a FACE which is normal to direction along given \a Axis;
 608 - an EDGE going along given \a Axis.
 609  \a i, \a j and \a k counts from zero.
 610
 611 Exception for \a Axis out of range.
 612 For 2D grids, \a k is a dummy argument. */
 613
 614 //================================================================================
 615 /*! Edge position to number conversion method*/
 616 //================================================================================
 617
 618 int GRID::getEdgeNumber(const int Axis, const int i, const int j, const int k)
 619   const throw (MEDEXCEPTION)
 620 {
 621   const char * LOC = "GRID::getEdgeNumber(Axis, i,j,k) :";
 622
 623   BEGIN_OF_MED(LOC);
 624
 625   int Len[4] = {0,_iArrayLength, _jArrayLength, _kArrayLength }, I=1, J=2, K=3;
 626   int maxAxis = Len[ K ] ? 3 : 2;
 627
 628   if (Axis <= 0 || Axis > maxAxis)
 629     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Axis out of range: " << Axis));
 630
 631   Len[Axis]--;
 632   int Nb = 1 + i + j*Len[ I ] + k*Len[ J ]*Len[ K ];
 633   Len[Axis]++ ;
 634
 635   if (!Len[ K ])
 636     Len[ K ] = 1;
 637
 638   if (Axis > 1) { // add all edges in i direction
 639     Len[I]-- ;
 640     Nb += Len[ I ]*Len[ J ]*Len[ K ];
 641     Len[I]++ ;
 642   }
 643   if (Axis > 2) { // add all edges in j direction
 644     Len[J]-- ;
 645     Nb += Len[ I ]*Len[ J ]*Len[ K ];
 646   }
 647
 648   END_OF_MED(LOC);
 649
 650   return Nb;
 651 }
 652
 653 //================================================================================
 654 /*!
 655 Returns a NODE, EDGE, FACE, CELL number by its position in the grid.
 656 Axis [1,2,3] means one of directions: along i, j or k
 657 For Cell contituents (FACE or EDGE), Axis selects one of those having same (i,j,k):
 658 - a FACE which is normal to direction along given Axis;
 659 - an EDGE going along given Axis.
 660 Exception for Axis out of range
 661 */
 662 //================================================================================
 663
 664 int GRID::getFaceNumber(const int Axis, const int i, const int j, const int k)
 665   const throw (MEDEXCEPTION)
 666 {
 667   const char * LOC = "GRID::getFaceNumber(Axis, i,j,k) :";
 668
 669   BEGIN_OF_MED(LOC);
 670
 671 //  if (Axis <= 0 || Axis > 3)
 672   if (Axis < 0 || Axis > 3)
 673     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Axis = " << Axis));
 674
 675   int Len[4] = {0,_iArrayLength-1, _jArrayLength-1, _kArrayLength-1 }, I=1, J=2, K=3;
 676
 677   Len[Axis]++;
 678   int Nb = 1 + i + j*Len[ I ] + k*Len[ I ]*Len[ J ];
 679   Len[Axis]--;
 680
 681   if (Axis > 1) // add all faces nornal to i direction
 682     Nb += ( Len[ I ]+1 )*Len[ J ]*Len[ K ];
 683
 684   if (Axis > 2) // add all faces nornal to j direction
 685     Nb += Len[ I ]*( Len[ J ]+1 )*Len[ K ];
 686
 687   END_OF_MED(LOC);
 688
 689   return Nb;
 690 }
 691 /*! @} */
 692
 693
 694 /*!
 695 @{
 696 \name Number to position conversion methods
 697
 698 \a getXXXPosition functions enable the user to convert
 699 a number into a \f$ (i,j,k) \f$ position.
 700            Axis [1,2,3] means one of directions: along i, j or k
 701            For Cell constituents (FACE or EDGE), Axis selects one of those having same (i,j,k):
 702            - a FACE which is normal to direction along given Axis;
 703            - an EDGE going along given Axis.
 704
 705     Exception for Number out of range.
 706 */
 707
 708 //================================================================================
 709 /*! Node number to position conversion method */
 710 //================================================================================
 711
 712 void GRID::getNodePosition(const int Number, int& i, int& j, int& k) const
 713   throw (MEDEXCEPTION)
 714 {
 715   const char * LOC = "GRID::getNodePosition(Number, i,j,k) :";
 716
 717   BEGIN_OF_MED(LOC);
 718
 719   if (Number <= 0 || Number > getNumberOfNodes() )
 720     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
 721
 722   int Len[] = {_iArrayLength, _jArrayLength, _kArrayLength }, I=0, J=1;
 723   // , K=2; !! UNUSED VARIABLE !!
 724
 725   int ijLen = Len[I] * Len[J]; // nb in a full k layer
 726   int kLen = (Number - 1) % ijLen;    // nb in the non full k layer
 727
 728   i = kLen % Len[J];
 729   j = kLen / Len[J];
 730   k = (Number - 1) / ijLen;
 731
 732   ////cout <<" NODE POS: " << Number << " - " << i << ", " << j << ", " << k << endl;
 733
 734   END_OF_MED(LOC);
 735
 736 }
 737
 738 //=======================================================================
 739 //function : getCellPosition
 740 //purpose  :
 741 //=======================================================================
 742 /*! Cell number to position conversion method */
 743 void GRID::getCellPosition(const int Number, int& i, int& j, int& k) const
 744   throw (MEDEXCEPTION)
 745 {
 746
 747   const char* LOC = "GRID::getCellPosition(Number, i,j,k) :";
 748   BEGIN_OF_MED(LOC);
 749
 750   int Len[4] = {0,_iArrayLength-1, _jArrayLength-1, _kArrayLength-1 }, I=1, J=2;
 751   // , K=3; !! UNUSED VARIABLE !!
 752
 753 //  if (Number <= 0 || Number > getCellNumber(Len[I]-1, Len[J]-1, Len[K]-1))
 754 //    throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
 755
 756   int ijLen = Len[I] * Len[J]; // nb in a full k layer
 757   int kLen = (Number - 1) % ijLen;    // nb in the non full k layer
 758
 759   i = kLen % Len[J];
 760   j = kLen / Len[J];
 761   k = (Number - 1) / ijLen;
 762
 763   END_OF_MED(LOC);
 764 }
 765
 766 //=======================================================================
 767 //function : getEdgePosition
 768 //purpose  :
 769 //=======================================================================
 770 /*! Edge number to poistion conversion method*/
 771 void GRID::getEdgePosition(const int Number, int& Axis, int& i, int& j, int& k)
 772   const throw (MEDEXCEPTION)
 773 {
 774   const char * LOC = "GRID::getEdgePosition(Number, i,j,k) :";
 775
 776   BEGIN_OF_MED(LOC);
 777
 778   if (!_jArrayLength)
 779     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "no edges in the grid: "));
 780
 781   if (Number <= 0)
 782     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
 783
 784   int Len[4] = {0,_iArrayLength, _jArrayLength, _kArrayLength }, I=1, J=2, K=3;
 785
 786   int theNb = Number;
 787   int maxAxis = _kArrayLength ? 3 : 2;
 788
 789   for (Axis = 1; Axis <= maxAxis; ++Axis)
 790   {
 791     Len[Axis]--;  // one less edge in Axis direction
 792
 793     // max edge number in Axis direction
 794     int maxNb = getEdgeNumber (Axis, Len[I]-1, Len[J]-1, Len[K]-1);
 795
 796     if (theNb > maxNb)
 797     {
 798       Len[Axis]++;
 799       theNb -= maxNb;
 800       continue;
 801     }
 802
 803     if (theNb == maxNb)
 804     {
 805       i = Len[I]-1;
 806       j = Len[J]-1;
 807       k = Len[K]-1;
 808     }
 809     else
 810     {
 811       int ijLen = Len[I] * Len[J]; // nb in a full k layer
 812       int kLen = (theNb - 1) % ijLen;    // nb in the non full k layer
 813       i = kLen % Len[J];
 814       j = kLen / Len[J];
 815       k = (theNb - 1) / ijLen;
 816     }
 817
 818   END_OF_MED(LOC);
 819
 820     return;
 821   }
 822
 823   throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
 824 }
 825
 826 //=======================================================================
 827 //function : getFacePosition
 828 //purpose  : return position (i,j,k) of an entity #Number
 829 //           Axis [1,2,3] means one of directions: along i, j or k
 830 //           For Cell contituents (FACE or EDGE), Axis selects one of those having same (i,j,k):
 831 //           * a FACE which is normal to direction along given Axis;
 832 //           * an EDGE going along given Axis.
 833 //           Exception for Number out of range
 834 //=======================================================================
 835 /*! Face number to position convertion method*/
 836 void GRID::getFacePosition(const int Number, int& Axis, int& i, int& j, int& k)
 837   const throw (MEDEXCEPTION)
 838 {
 839   const char * LOC = "GRID::getFacePosition(Number, i,j,k) :";
 840
 841   BEGIN_OF_MED(LOC);
 842
 843   if (_kArrayLength == 0) {
 844     getCellPosition(Number, i, j, k);
 845     Axis = 1;
 846     return;
 847   };
 848
 849   if (!_kArrayLength)
 850     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "no faces in the grid: "));
 851
 852   if ( Number <= 0 )
 853     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
 854
 855   int Len[4] = {0,_iArrayLength-1, _jArrayLength-1, _kArrayLength-1 }, I=1, J=2, K=3;
 856   int theNb = Number;
 857
 858   for (Axis = 1; Axis <= 3; ++Axis)
 859   {
 860     Len[Axis]++;
 861
 862     // max face number in Axis direction
 863     int maxNb = getFaceNumber (Axis, Len[I]-1, Len[J]-1, Len[K]-1);
 864
 865     if (theNb > maxNb)
 866     {
 867       Len[Axis]--;
 868       theNb -= maxNb;
 869       continue;
 870     }
 871
 872     if (theNb == maxNb)
 873     {
 874       i = Len[I]-1;
 875       j = Len[J]-1;
 876       k = Len[K]-1;
 877     }
 878     else
 879     {
 880       int ijLen = Len[I] * Len[J]; // nb in a full k layer
 881       int kLen = (theNb - 1)  % ijLen;    // nb in the non full k layer
 882       i = kLen % Len[J];
 883       j = kLen / Len[J];
 884       k = (theNb - 1) / ijLen;
 885     }
 886
 887   END_OF_MED(LOC);
 888
 889     return;
 890   }
 891
 892   throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
 893 }
 894 /*!
 895 @}
 896 \if MEDMEM_ug
 897 @}
 898 \endif
 899 */
 900
 901 //================================================================================
 902 /*! Get the number of different geometric types for a given entity type.
 903
 904     medEntityMesh entity : MED_CELL, MED_FACE, MED_EDGE, MED_NODE,
 905     MED_ALL_ENTITIES
 906
 907 */
 908 //================================================================================
 909
 910 int GRID::getNumberOfTypes(MED_EN::medEntityMesh entity) const
 911 {
 912   MESSAGE_MED("GRID::getNumberOfTypes(medEntityMesh entity) : "<<entity);
 913   return 1; // a grid has one type
 914 }
 915
 916 //================================================================================
 917 /*!
 918   Return the number of element of given geometric type of given entity. Return 0 if query is not defined.
 919 */
 920 //================================================================================
 921
 922 int GRID::getNumberOfElements(MED_EN::medEntityMesh entity, MED_EN::medGeometryElement Type) const
 923 {
 924   int numberOfElements=0;
 925
 926   // Cas où le nombre d'éléments n'est pas nul
 927   if (entity==MED_EN::MED_FACE && (Type==MED_EN::MED_QUAD4 || Type==MED_EN::MED_ALL_ELEMENTS) && getMeshDimension()>2)
 928     numberOfElements=
 929       (_iArrayLength-1)*(_jArrayLength-1)*(_kArrayLength  )+
 930       (_iArrayLength-1)*(_jArrayLength  )*(_kArrayLength-1)+
 931       (_iArrayLength  )*(_jArrayLength-1)*(_kArrayLength-1);
 932
 933   else if (entity==MED_EN::MED_EDGE && (Type==MED_EN::MED_SEG2 || Type==MED_EN::MED_ALL_ELEMENTS))
 934     if ( _spaceDimension==2)
 935       numberOfElements=_iArrayLength*(_jArrayLength-1) + (_iArrayLength-1)*_jArrayLength;
 936     else if ( _spaceDimension==1)
 937       numberOfElements=_iArrayLength-1;
 938     else // 3D
 939       numberOfElements=
 940         (_iArrayLength*(_jArrayLength-1) + (_iArrayLength-1)*_jArrayLength) * _kArrayLength +
 941         _iArrayLength*_jArrayLength*(_kArrayLength-1);
 942
 943   else if (entity==MED_EN::MED_NODE && (Type==MED_EN::MED_NONE || Type==MED_EN::MED_ALL_ELEMENTS) && _spaceDimension>0)
 944     numberOfElements=getNumberOfNodes();
 945
 946   else if (entity==MED_EN::MED_CELL && _spaceDimension==3 && (Type==MED_EN::MED_HEXA8 || Type==MED_EN::MED_ALL_ELEMENTS) )
 947     numberOfElements=(_iArrayLength-1)*(_jArrayLength-1)*(_kArrayLength-1);
 948
 949   else if (entity==MED_EN::MED_CELL && _spaceDimension==2 && (Type==MED_EN::MED_QUAD4 || Type==MED_EN::MED_ALL_ELEMENTS))
 950     numberOfElements=(_iArrayLength-1)*(_jArrayLength-1);
 951
 952   else if (entity==MED_EN::MED_CELL && _spaceDimension==1 && (Type==MED_EN::MED_SEG2 || Type==MED_EN::MED_ALL_ELEMENTS) )
 953     numberOfElements=_iArrayLength-1;
 954
 955   MESSAGE_MED("GRID::getNumberOfElements - entity=" << entity << " Type=" << Type);
 956   MESSAGE_MED("_spaceDimension=" << _spaceDimension << "  numberOfElements=" << numberOfElements);
 957
 958   return numberOfElements;
 959 }
 960
 961 //================================================================================
 962 /*!
 963   Return the geometric type of global element Number of entity Entity.
 964 */
 965 //================================================================================
 966
 967 MED_EN::medGeometryElement GRID::getElementType(MED_EN::medEntityMesh Entity,int Number) const
 968 {
 969   return getTypes(Entity)[0];
 970 }
 971
 972 //================================================================================
 973 /*!
 974  * \brief Return mesh dimension
 975  */
 976 //================================================================================
 977
 978 int GRID::getMeshDimension() const
 979 {
 980   return getSpaceDimension();
 981 }
 982
 983 //================================================================================
 984 /*!
 985  * \brief It is a grid
 986  */
 987 //================================================================================
 988
 989 bool GRID::getIsAGrid() const
 990 {
 991   return true;
 992 }
 993
 994 //================================================================================
 995 /*!
 996  * \brief Return number of nodes
 997  */
 998 //================================================================================
 999
1000 int GRID::getNumberOfNodes() const
1001 {
1002   if ( _gridType == MED_EN::MED_BODY_FITTED )
1003     return _coordinate ? _coordinate->getNumberOfNodes() : 0;
1004
1005   switch ( _spaceDimension )
1006   {
1007   case 3: return _iArrayLength * _jArrayLength * _kArrayLength;
1008   case 2: return _iArrayLength * _jArrayLength;
1009   case 1: return _iArrayLength;
1010   }
1011   return 0;
1012 }
1013
1014 //=======================================================================
1015 /*!
1016  * Returns "CARTESIAN", "CYLINDRICAL" or "SPHERICAL"
1017  */
1018 //=======================================================================
1019
1020 std::string         GRID::getCoordinatesSystem() const
1021 {
1022   return _coordinate ? _coordinate->getCoordinatesSystem() : defaultStrings()[0];
1023 }
1024
1025 //=======================================================================
1026 /*!
1027  * Returns an array with names of coordinates. \n
1028  *     Example : \n
1029  *     - x,y,z \n
1030  *     - r,teta,phi \n
1031  *     - ...
1032  */
1033 //=======================================================================
1034
1035 const std::string * GRID::getCoordinatesNames() const
1036 {
1037   return _coordinate ? _coordinate->getCoordinatesNames() : defaultStrings();
1038 }
1039
1040 //=======================================================================
1041 /*!
1042  * Returns an array with units of coordinates (cm, m, mm, ...)
1043  * It could be empty. We suppose we are IS (meter).
1044  */
1045 //=======================================================================
1046
1047 const std::string * GRID::getCoordinatesUnits() const
1048 {
1049   return _coordinate ? _coordinate->getCoordinatesUnits() : defaultStrings();
1050 }
1051
1052 //=======================================================================
1053 /*!
1054   Returns a support which reference all elements on the boundary of mesh.
1055   For a d-dimensional mesh, a boundary element is defined as a d-1 dimension
1056   element that is referenced by only one element in the full descending connectivity.
1057
1058   This method can also return the list of nodes that belong to the boundary elements.
1059
1060   WARNING: This method can recalculate descending connectivity from partial to full form,
1061   so that partial SUPPORT on d-1 dimension elements becomes invalid.
1062
1063   \param Entity entity on which the boundary is desired. It has to be either \a MED_NODE or the
1064   d-1 dimension entity type (MED_FACE in 3D, MED_EDGE in 2D).
1065 */
1066 //=======================================================================
1067
1068 SUPPORT * GRID::getBoundaryElements(MED_EN::medEntityMesh Entity) const throw (MEDEXCEPTION)
1069 {
1070   const char * LOC = "GMESH::getBoundaryElements() : " ;
1071   if ( _spaceDimension < 2 )
1072     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not implemented in " << _spaceDimension <<"D space !"));
1073
1074   if ( _gridType == MED_POLAR)
1075     throw MEDEXCEPTION("GRID::getBoundaryElements() : not implemented on MED_POLAR grig");
1076
1077
1078   medEntityMesh entityToParse=Entity;
1079   if(_spaceDimension == 3)
1080     if (Entity != MED_FACE)
1081       {
1082         if(Entity==MED_NODE)
1083           entityToParse=MED_FACE;
1084         else
1085           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<Entity<<" !"));
1086       }
1087   if(_spaceDimension == 2)
1088     if(Entity != MED_EDGE)
1089       {
1090         if(Entity==MED_NODE)
1091           entityToParse=MED_EDGE;
1092         else
1093           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<Entity<<" !"));
1094       }
1095   list<int> bnd_elems1, bnd_elems2;
1096   int numberOf = getNumberOfElements(entityToParse,MED_ALL_ELEMENTS) ;
1097
1098   if ( _coordinate->getNumberOfNodes() > 0 ) // BODY FITTED
1099     {
1100       throw MEDEXCEPTION("GRID::getBoundaryElements() : not implemented on BOBY FITTED grig");
1101     }
1102   else if ( entityToParse == MED_FACE ) // 3D CARTESIAN
1103     {
1104       const int nb_x = getArrayLength( 1 ) - 1;
1105       const int nb_y = getArrayLength( 2 ) - 1;
1106       const int nb_z = getArrayLength( 3 ) - 1;
1107       // normal to OX
1108       for ( int z = 0; z < nb_z; ++z )
1109         for ( int y = 0; y < nb_y; ++y )
1110           {
1111             bnd_elems1.push_back( getFaceNumber( 1, 0,    y, z ));
1112             bnd_elems2.push_back( getFaceNumber( 1, nb_x, y, z ));
1113           }
1114       bnd_elems1.splice( bnd_elems1.end(), bnd_elems2 ); // to have ids in increasing order
1115
1116       // normal to OY
1117       for ( int z = 0; z < nb_z; ++z )
1118         for ( int x = 0; x < nb_x; ++x )
1119           {
1120             bnd_elems1.push_back( getFaceNumber( 2, x, 0,    z ));
1121             bnd_elems2.push_back( getFaceNumber( 2, x, nb_y, z ));
1122           }
1123       bnd_elems1.splice( bnd_elems1.end(), bnd_elems2 );
1124
1125       // normal to OZ
1126       for ( int y = 0; y < nb_y; ++y )
1127         for ( int x = 0; x < nb_x; ++x )
1128           {
1129             bnd_elems1.push_back( getFaceNumber( 3, x, y, 0    ));
1130             bnd_elems2.push_back( getFaceNumber( 3, x, y, nb_z ));
1131           }
1132       bnd_elems1.splice( bnd_elems1.end(), bnd_elems2 );
1133     }
1134   else
1135     {
1136       const int nb_x = getArrayLength( 1 ) - 1;
1137       const int nb_y = getArrayLength( 2 ) - 1;
1138       // edge || OX
1139       for ( int x = 0; x < nb_x; ++x )
1140         {
1141           bnd_elems1.push_back( getEdgeNumber( 1, x, 0 ));
1142           bnd_elems2.push_back( getEdgeNumber( 1, x, nb_y ));
1143         }
1144       bnd_elems1.splice( bnd_elems1.end(), bnd_elems2 ); // to have ids in increasing order
1145       // edge || OY
1146       for ( int y = 0; y < nb_y; ++y )
1147         {
1148           bnd_elems1.push_back( getEdgeNumber( 2, y, 0 ));
1149           bnd_elems2.push_back( getEdgeNumber( 2, y, nb_x ));
1150         }
1151       bnd_elems1.splice( bnd_elems1.end(), bnd_elems2 );
1152     }
1153
1154   if ( bnd_elems1.empty() && numberOf != 0 )
1155     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No boundary elements found by reverse descending connectivity for entity "<<Entity<<" !"));
1156
1157   if ( Entity == MED_NODE )
1158     return buildSupportOnNodeFromElementList(bnd_elems1,entityToParse);
1159   else
1160     return buildSupportOnElementsFromElementList(bnd_elems1,entityToParse);
1161 }
1162
1163 SUPPORT * GRID::getSkin(const SUPPORT * Support3D) throw (MEDEXCEPTION)
1164 {
1165   throw MEDEXCEPTION("GRID::getSkin() : Not implemented yet");
1166 }
1167 SUPPORT *GRID::buildSupportOnNodeFromElementList(const std::list<int>& listOfElt,
1168                                                  MED_EN::medEntityMesh entity) const
1169   throw (MEDEXCEPTION)
1170 {
1171   throw MEDEXCEPTION("GRID::buildSupportOnNodeFromElementList() : Not implemented yet");
1172 }
1173 void GRID::fillSupportOnNodeFromElementList(const std::list<int>& listOfElt,
1174                                             SUPPORT *             supportToFill) const
1175   throw (MEDEXCEPTION)
1176 {
1177   throw MEDEXCEPTION("GRID::fillSupportOnNodeFromElementList() : Not implemented yet");
1178 }
1179
1180 FIELD<double>* GRID::getVolume (const SUPPORT * Support, bool isAbs) const
1181   throw (MEDEXCEPTION)
1182 {
1183   throw MEDEXCEPTION("GRID::getVolume() : Not implemented yet");
1184 }
1185
1186 FIELD<double>* GRID::getArea (const SUPPORT * Support) const
1187   throw (MEDEXCEPTION)
1188 {
1189   throw MEDEXCEPTION("GRID::getArea() : Not implemented yet");
1190 }
1191
1192 FIELD<double>* GRID::getLength (const SUPPORT * Support) const
1193   throw (MEDEXCEPTION)
1194 {
1195   throw MEDEXCEPTION("GRID::getLength() : Not implemented yet");
1196 }
1197
1198 FIELD<double>* GRID::getNormal (const SUPPORT * Support) const
1199   throw (MEDEXCEPTION)
1200 {
1201   throw MEDEXCEPTION("GRID::getNormal() : Not implemented yet");
1202 }
1203
1204 FIELD<double>* GRID::getBarycenter (const SUPPORT * Support) const
1205   throw (MEDEXCEPTION)
1206 {
1207   throw MEDEXCEPTION("GRID::getBarycenter() : Not implemented yet");
1208 }
1209
1210 vector< vector<double> >   GRID::getBoundingBox() const
1211 {
1212   throw MEDEXCEPTION("GRID::getBoundingBox() : Not implemented yet");
1213 }