Back to index

salome-med  6.5.0
MEDMEM_Grid.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00004 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 
00023 // File      : MEDMEM_Grid.hxx
00024 // Created   : Wed Dec 18 08:35:26 2002
00025 // Descr     : class containing structured mesh data
00026 // Author    : Edward AGAPOV (eap)
00027 // Project   : SALOME Pro
00028 // Module    : MED 
00029 //
00030 #include "MEDMEM_Grid.hxx"
00031 #include "MEDMEM_Meshing.hxx"
00032 #include "MEDMEM_CellModel.hxx"
00033 #include "MEDMEM_SkyLineArray.hxx"
00034 #include "MEDMEM_DriverFactory.hxx"
00035 
00036 #include <memory>
00037 
00038 using namespace std;
00039 using namespace MEDMEM;
00040 using namespace MED_EN;
00041 
00042 
00043 // Block defining the comments for the MEDMEM_ug documentation
00044 
00058 namespace
00059 {
00060   const string* defaultStrings()
00061   {
00062     static const string str[3] = { "UNDEFINED", "UNDEFINED", "UNDEFINED" };
00063     return str;
00064   }
00065 }
00066 
00067 //=======================================================================
00068 //function : GRID
00069 //purpose  : empty constructor
00070 //=======================================================================
00071 
00072 GRID::GRID() {
00073   init();
00074   MESSAGE_MED("A GRID CREATED");
00075 }
00076 //
00077 //=======================================================================
00081 
00099 GRID::GRID(const std::vector<std::vector<double> >& xyz_array,
00100            const std::vector<std::string>&          coord_name,
00101            const std::vector<std::string>&          coord_unit,
00102            const MED_EN::med_grid_type              type)
00103   :_gridType(type)
00104 {
00105     init(); // PAL 12136
00106     _is_default_gridType = false;
00107 
00108     _spaceDimension = xyz_array.size();
00109 
00110     // create a non allocated COORDINATE
00111     _coordinate = new COORDINATE(_spaceDimension, &coord_name[0], &coord_unit[0]);
00112     string coordinateSystem = "UNDEFINED";
00113     if( _gridType == MED_CARTESIAN) 
00114       coordinateSystem = "CARTESIAN";
00115     else if ( _gridType == MED_POLAR) 
00116       coordinateSystem = "CYLINDRICAL";
00117     _coordinate->setCoordinatesSystem(coordinateSystem);
00118 
00119     // set the GRID part
00120     if (_spaceDimension>=1)
00121     {
00122         _iArrayLength=xyz_array[0].size();
00123         _iArray=new double[_iArrayLength];
00124         std::copy(xyz_array[0].begin(),xyz_array[0].end(),_iArray);
00125     }
00126     if (_spaceDimension>=2)
00127     {
00128         _jArrayLength=xyz_array[1].size();
00129         _jArray=new double[_jArrayLength];
00130         std::copy(xyz_array[1].begin(),xyz_array[1].end(),_jArray);
00131     }
00132     if (_spaceDimension>=3)
00133     {
00134         _kArrayLength=xyz_array[2].size();
00135         _kArray=new double[_kArrayLength];
00136         std::copy(xyz_array[2].begin(),xyz_array[2].end(),_kArray);
00137     }
00138 }
00139 
00140 //================================================================================
00147 //================================================================================
00148 
00149 GRID::GRID(driverTypes driverType, const string &  fileName, const string &  driverName)
00150 {
00151   
00152   const char* LOC = "GRID::GRID(driverTypes , const string &  , const string &) : ";
00153   BEGIN_OF_MED(LOC);
00154   
00155   init();
00156   auto_ptr<GENDRIVER> myDriver (DRIVERFACTORY::buildDriverForMesh(driverType,fileName,this,driverName,RDONLY));
00157   myDriver->open();
00158   myDriver->read();
00159   myDriver->close();
00160 
00161   END_OF_MED(LOC);
00162 }
00163 
00166 //=======================================================================
00167 //function : GRID
00168 //purpose  : empty constructor
00169 //=======================================================================
00170 
00171 GRID::GRID(const MED_EN::med_grid_type type)
00172 {
00173   init();
00174   _gridType = type;
00175   _is_default_gridType = false;
00176   MESSAGE_MED("A TYPED GRID CREATED");
00177 }
00178 
00179 //=======================================================================
00180 //function : GRID
00181 //purpose  : copying constructor
00182 //=======================================================================
00183 
00184 GRID::GRID(const GRID& otherGrid) {
00185   *this = otherGrid;
00186 }
00187 
00188 //=======================================================================
00189 //function : ~GRID
00190 //purpose  : destructor
00191 //=======================================================================
00192 
00193 GRID::~GRID() {
00194   MESSAGE_MED("GRID::~GRID() : Destroying the Grid");
00195   if ( _coordinate ) delete _coordinate; _coordinate = 0;
00196   if ( _iArray != (double* ) NULL) delete [] _iArray;
00197   if ( _jArray != (double* ) NULL) delete [] _jArray;
00198   if ( _kArray != (double* ) NULL) delete [] _kArray;
00199 }
00200 
00201 //=======================================================================
00202 //function : init
00203 //purpose : 
00204 //=======================================================================
00205 
00206 void GRID::init()
00207 {
00208   GMESH::init();
00209 
00210   _gridType = MED_CARTESIAN;
00211   _is_default_gridType = true;
00212   _coordinate = 0;
00213   _iArray = _jArray = _kArray = (double* ) NULL;
00214   _iArrayLength = _jArrayLength = _kArrayLength = 0;
00215 
00216 }
00217 
00218 //================================================================================
00222 //================================================================================
00223 
00224 bool GRID::isEmpty() const
00225 {
00226   return !_iArrayLength && !_coordinate;
00227 }
00228 
00229 
00230 //=======================================================================
00231 //function: operator=
00232 //purpose : operator=
00233 //=======================================================================
00234 
00235 GRID & GRID::operator=(const GRID & otherGrid)
00236 {
00237   // NOT IMPLEMENTED
00238 
00239   GMESH::operator=(otherGrid);
00240   return *this;
00241 }
00242 
00243 //=======================================================================
00250 //=======================================================================
00251 
00252 bool GRID::deepCompare(const GMESH& gother) const
00253 {
00254   if ( gother.getIsAGrid() != getIsAGrid())
00255     return false;
00256   const GRID& other = static_cast<const GRID&>( gother );
00257 
00258   if ( getSpaceDimension() != other.getSpaceDimension() )
00259     return false;
00260 
00261   if ( _gridType != other._gridType )
00262     return false;
00263 
00264   if( bool( _coordinate) != bool(other._coordinate))
00265     return false;
00266 
00267   if ( _coordinate->getNumberOfNodes() > 0 )
00268   {
00269     if ( _coordinate->getNumberOfNodes() != other._coordinate->getNumberOfNodes() )
00270       return false;
00271     const int size = _coordinate->getNumberOfNodes() * getSpaceDimension();
00272     const double* coord1=_coordinate->getCoordinates(MED_FULL_INTERLACE);
00273     const double* coord2=other._coordinate->getCoordinates(MED_FULL_INTERLACE);
00274     for(int i = 0 ; i < size; i++ )
00275       if ( !(fabs(coord1[i]-coord2[i])<1e-15))
00276         return false;
00277   }
00278 
00279   if ( _iArrayLength != other._iArrayLength )
00280     return false;
00281   if ( _jArrayLength != other._jArrayLength )
00282     return false;
00283   if ( _kArrayLength != other._kArrayLength )
00284     return false;
00285 
00286   if ( bool(_iArray) != bool(other._iArray) )
00287     return false;
00288   if ( bool(_jArray) != bool(other._jArray) )
00289     return false;
00290   if ( bool(_kArray) != bool(other._kArray) )
00291     return false;
00292 
00293   if ( _iArray )
00294     for ( int i = 0; i < _iArrayLength; ++i )
00295       if ( !(fabs(_iArray[i]-other._iArray[i])<1e-15))
00296         return false;
00297   if ( _jArray )
00298     for ( int i = 0; i < _jArrayLength; ++i )
00299       if ( !(fabs(_jArray[i]-other._jArray[i])<1e-15))
00300         return false;
00301   if ( _kArray )
00302     for ( int i = 0; i < _kArrayLength; ++i )
00303       if ( !(fabs(_kArray[i]-other._kArray[i])<1e-15))
00304         return false;
00305 
00306   return true;
00307 }
00308 
00309 //================================================================================
00313 //================================================================================
00314 
00315 void GRID::printMySelf(std::ostream &os) const
00316 {
00317   // TODO
00318   cout << "NOT IMPLEMENTED" << endl;
00319 }
00320 
00321 //================================================================================
00325 //================================================================================
00326 
00327 const MESH * GRID::convertInMESH() const
00328 {
00329   MESHING* mesh = new MESHING;
00330   mesh->setName( getName() );
00331 
00332   int i, j, k;
00333 
00334   // ---------------
00335   // 1. Coordinates
00336   // ---------------
00337 
00338   PointerOf< double > coords;
00339   if ( _gridType == MED_BODY_FITTED )
00340   {
00341     if ( !_coordinate )
00342       throw MEDEXCEPTION
00343         (LOCALIZED("GRID::convertInMESH() : _coordinate of MED_BODY_FITTED GRID not defined !"));
00344     coords.set( _coordinate->getCoordinates(MED_FULL_INTERLACE));
00345   }
00346   else
00347   {
00348     coords.set( getSpaceDimension() * getNumberOfNodes() );
00349     double* coord = coords;
00350     switch ( getSpaceDimension() )
00351     {
00352     case 3:
00353       for (k=0; k < _kArrayLength; ++k) 
00354         for (j=0; j < _jArrayLength; ++j)
00355           for (i=0; i < _iArrayLength; ++i)
00356             *coord++ = _iArray[i], *coord++ = _jArray[j], *coord++ = _kArray[k];
00357       break;
00358 
00359     case 2:
00360       for (j=0; j < _jArrayLength; ++j)
00361         for (i=0; i < _iArrayLength; ++i)
00362           *coord++ = _iArray[i], *coord++ = _jArray[j];
00363       break;
00364 
00365     case 1:
00366       coords.set(_iArray);
00367       break;
00368     }
00369   }
00370   mesh->setCoordinates( getSpaceDimension(),
00371                         getNumberOfNodes(),
00372                         (const double *) coords,
00373                         getCoordinatesSystem(),
00374                         MED_EN::MED_FULL_INTERLACE);
00375   mesh->setCoordinatesNames( getCoordinatesNames() );
00376   mesh->setCoordinatesUnits( getCoordinatesUnits() );
00377 
00378   // ----------------
00379   // 2. Connectivity
00380   // ----------------
00381 
00382   // 2.1 Cells
00383   // ---------
00384   medEntityMesh subEntity;
00385   {
00386     mesh->setNumberOfTypes( getNumberOfTypes(MED_CELL), MED_CELL );
00387     mesh->setTypes( getTypes( MED_CELL), MED_CELL );
00388     int nbCells = getNumberOfElements( MED_CELL, MED_ALL_ELEMENTS );
00389     mesh->setNumberOfElements( &nbCells, MED_CELL );
00390 
00391     vector<int> conn;
00392     switch ( _spaceDimension )
00393     {
00394     case 3: // HEXA8
00395       for ( k = 0; k < _kArrayLength-1; ++k )
00396         for ( j = 0; j < _jArrayLength-1; ++j )
00397           for ( i = 0; i < _iArrayLength-1; ++i )
00398           {
00399             conn.push_back( getNodeNumber( i  , j  , k  ));
00400             conn.push_back( getNodeNumber( i  , j+1, k  ));
00401             conn.push_back( getNodeNumber( i+1, j+1, k  ));
00402             conn.push_back( getNodeNumber( i+1, j  , k  ));
00403             conn.push_back( getNodeNumber( i  , j  , k+1));
00404             conn.push_back( getNodeNumber( i  , j+1, k+1));
00405             conn.push_back( getNodeNumber( i+1, j+1, k+1));
00406             conn.push_back( getNodeNumber( i+1, j  , k+1));
00407           }
00408       subEntity = MED_FACE;
00409       break;
00410 
00411     case 2: // QUAD4
00412       for ( j = 0; j < _jArrayLength-1; ++j )
00413         for ( i = 0; i < _iArrayLength-1; ++i )
00414         {
00415           int n1 = 1 + i + j*_iArrayLength;
00416           conn.push_back( n1 );
00417           conn.push_back( n1 + _iArrayLength );
00418           conn.push_back( n1 + _iArrayLength + 1 );
00419           conn.push_back( n1 + 1 );
00420         }
00421       subEntity = MED_EDGE;
00422       break;
00423 
00424     case 1: // SEG2
00425       for ( i = 0; i < _iArrayLength-1; ++i )
00426       {
00427         conn.push_back( i + 1 );
00428         conn.push_back( i + 2 );
00429       }
00430       break;
00431     }
00432     mesh->setConnectivity( MED_CELL, getTypes(MED_CELL)[0], &conn[0] );
00433   }
00434 
00435   // 2.2 sub entities
00436   // -----------------
00437 
00438   if ( _spaceDimension > 1 )
00439   {
00440     mesh->setNumberOfTypes( getNumberOfTypes(subEntity), subEntity );
00441     mesh->setTypes( getTypes( subEntity), subEntity );
00442     int nbCells = getNumberOfElements( subEntity, MED_ALL_ELEMENTS );
00443     mesh->setNumberOfElements( &nbCells, subEntity );
00444 
00445     vector<int> conn;
00446     switch ( _spaceDimension )
00447     {
00448     case 3: // QUAD4
00449 
00450       // normal to OX
00451       for ( k = 0; k < _kArrayLength-1; ++k )
00452         for ( j = 0; j < _jArrayLength-1; ++j )
00453           for ( i = 0; i < _iArrayLength; ++i )
00454           {
00455             conn.push_back( getNodeNumber( i, j  , k  ));
00456             conn.push_back( getNodeNumber( i, j  , k+1));
00457             conn.push_back( getNodeNumber( i, j+1, k+1));
00458             conn.push_back( getNodeNumber( i, j+1, k  ));
00459           }
00460       // normal to OY
00461       for ( k = 0; k < _kArrayLength-1; ++k )
00462         for ( j = 0; j < _jArrayLength; ++j )
00463           for ( i = 0; i < _iArrayLength-1; ++i )
00464           {
00465             conn.push_back( getNodeNumber( i  , j, k  ));
00466             conn.push_back( getNodeNumber( i+1, j, k  ));
00467             conn.push_back( getNodeNumber( i+1, j, k+1));
00468             conn.push_back( getNodeNumber( i  , j, k+1));
00469           }
00470       // normal to OZ
00471       for ( k = 0; k < _kArrayLength; ++k )
00472         for ( j = 0; j < _jArrayLength-1; ++j )
00473           for ( i = 0; i < _iArrayLength-1; ++i )
00474           {
00475             conn.push_back( getNodeNumber( i  , j  , k));
00476             conn.push_back( getNodeNumber( i  , j+1, k));
00477             conn.push_back( getNodeNumber( i+1, j+1, k));
00478             conn.push_back( getNodeNumber( i+1, j  , k));
00479           }
00480       break;
00481 
00482     case 2: // SEG2
00483 
00484       // || OX
00485       for ( j = 0; j < _jArrayLength; ++j )
00486         for ( i = 0; i < _iArrayLength-1; ++i )
00487         {
00488           int n1 = 1 + i + j*_iArrayLength;
00489           conn.push_back( n1 );
00490           conn.push_back( n1 + 1 );
00491         }
00492       // || OY
00493       for ( j = 0; j < _jArrayLength-1; ++j )
00494         for ( i = 0; i < _iArrayLength; ++i )
00495         {
00496           int n1 = 1 + i + j*_iArrayLength;
00497           conn.push_back( n1 );
00498           conn.push_back( n1 + _iArrayLength );
00499         }
00500       break;
00501 
00502     }
00503     mesh->setConnectivity( subEntity, getTypes(subEntity)[0], &conn[0] );
00504   }
00505 
00506   // 3. Groups and Families
00507   // -----------------------
00508 
00509   const vector<GROUP*>* groups[] = { &_groupNode, &_groupCell, &_groupFace, &_groupEdge };
00510   for ( i = 0; i < 4; ++i )
00511     for ( j = 0; j < (int)groups[i]->size(); ++j )
00512       mesh->addGroup( * groups[i]->at( j ));
00513 
00514   return mesh;
00515 }
00516 
00517 //=======================================================================
00521 //=======================================================================
00522 
00523 const medGeometryElement * GRID::getTypes(MED_EN::medEntityMesh entity) const
00524 {
00525   static const medGeometryElement _gridGeometry[4]={MED_HEXA8,MED_QUAD4,MED_SEG2,MED_POINT1};
00526   int i=0;
00527   if(entity==MED_CELL)
00528   {
00529     i=3-_spaceDimension;
00530   }
00531   else if(entity==MED_FACE && _spaceDimension>2 )
00532     i=1;
00533   else if(entity==MED_EDGE && _spaceDimension>1 )
00534     i=2;
00535   else if(entity==MED_NODE && _spaceDimension>0)
00536     i=3;
00537   else
00538     throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::getGeometricTypes : Entity not defined !"));
00539   return &_gridGeometry[i];
00540 }
00541 
00542 //=======================================================================
00543 //function : getArrayLength
00544 //purpose  : return array length. Axis = [1,2,3] meaning [i,j,k],
00545 //=======================================================================
00553 int GRID::getArrayLength( const int Axis ) const throw (MEDEXCEPTION)
00554 {
00555   switch (Axis) {
00556   case 1: return _iArrayLength;
00557   case 2: return _jArrayLength;
00558   case 3: return _kArrayLength;
00559   default:
00560     throw MED_EXCEPTION ( LOCALIZED( STRING("GRID::getArrayLength ( ") << Axis << ")"));
00561   }
00562   return 0;
00563 }
00564 
00565 //=======================================================================
00566 //function : getArrayValue
00567 //purpose  : return i-th array component. Axis = [1,2,3] meaning [i,j,k],
00568 //           exception if Axis out of [1-3] range
00569 //           exception if i is out of range 0 <= i < getArrayLength(Axis);
00570 //=======================================================================
00574 const double GRID::getArrayValue (const int Axis, const int i) const throw (MEDEXCEPTION)
00575 {
00576   if (i < 0 || i >= getArrayLength(Axis))
00577     throw MED_EXCEPTION
00578       ( LOCALIZED(STRING("GRID::getArrayValue ( ") << Axis << ", " << i << ")"));
00579   
00580   switch (Axis) {
00581   case 1:  return _iArray[ i ];
00582   case 2:  return _jArray[ i ];
00583   default: return _kArray[ i ];
00584   }
00585   return 0.0;
00586 }
00587 
00614 //================================================================================
00616 //================================================================================
00617 
00618 int GRID::getEdgeNumber(const int Axis, const int i, const int j, const int k)
00619   const throw (MEDEXCEPTION)
00620 {
00621   const char * LOC = "GRID::getEdgeNumber(Axis, i,j,k) :";
00622 
00623   BEGIN_OF_MED(LOC);
00624 
00625   int Len[4] = {0,_iArrayLength, _jArrayLength, _kArrayLength }, I=1, J=2, K=3;
00626   int maxAxis = Len[ K ] ? 3 : 2;
00627   
00628   if (Axis <= 0 || Axis > maxAxis)
00629     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Axis out of range: " << Axis));
00630 
00631   Len[Axis]--;
00632   int Nb = 1 + i + j*Len[ I ] + k*Len[ J ]*Len[ K ];
00633   Len[Axis]++ ;
00634   
00635   if (!Len[ K ])
00636     Len[ K ] = 1; 
00637 
00638   if (Axis > 1) { // add all edges in i direction
00639     Len[I]-- ;
00640     Nb += Len[ I ]*Len[ J ]*Len[ K ];
00641     Len[I]++ ;
00642   }
00643   if (Axis > 2) { // add all edges in j direction
00644     Len[J]-- ;
00645     Nb += Len[ I ]*Len[ J ]*Len[ K ];
00646   }
00647 
00648   END_OF_MED(LOC);
00649 
00650   return Nb;
00651 }
00652 
00653 //================================================================================
00662 //================================================================================
00663 
00664 int GRID::getFaceNumber(const int Axis, const int i, const int j, const int k)
00665   const throw (MEDEXCEPTION)
00666 {
00667   const char * LOC = "GRID::getFaceNumber(Axis, i,j,k) :";
00668   
00669   BEGIN_OF_MED(LOC);
00670 
00671 //  if (Axis <= 0 || Axis > 3)
00672   if (Axis < 0 || Axis > 3)
00673     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Axis = " << Axis));
00674 
00675   int Len[4] = {0,_iArrayLength-1, _jArrayLength-1, _kArrayLength-1 }, I=1, J=2, K=3;
00676 
00677   Len[Axis]++;
00678   int Nb = 1 + i + j*Len[ I ] + k*Len[ I ]*Len[ J ];
00679   Len[Axis]--;
00680 
00681   if (Axis > 1) // add all faces nornal to i direction
00682     Nb += ( Len[ I ]+1 )*Len[ J ]*Len[ K ];
00683 
00684   if (Axis > 2) // add all faces nornal to j direction
00685     Nb += Len[ I ]*( Len[ J ]+1 )*Len[ K ];
00686 
00687   END_OF_MED(LOC);
00688 
00689   return Nb;
00690 }
00708 //================================================================================
00710 //================================================================================
00711 
00712 void GRID::getNodePosition(const int Number, int& i, int& j, int& k) const
00713   throw (MEDEXCEPTION)
00714 {
00715   const char * LOC = "GRID::getNodePosition(Number, i,j,k) :";
00716   
00717   BEGIN_OF_MED(LOC);
00718 
00719   if (Number <= 0 || Number > getNumberOfNodes() )
00720     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
00721 
00722   int Len[] = {_iArrayLength, _jArrayLength, _kArrayLength }, I=0, J=1;
00723   // , K=2; !! UNUSED VARIABLE !!
00724 
00725   int ijLen = Len[I] * Len[J]; // nb in a full k layer
00726   int kLen = (Number - 1) % ijLen;    // nb in the non full k layer
00727   
00728   i = kLen % Len[J];
00729   j = kLen / Len[J];
00730   k = (Number - 1) / ijLen;
00731 
00733 
00734   END_OF_MED(LOC);
00735 
00736 }
00737 
00738 //=======================================================================
00739 //function : getCellPosition
00740 //purpose  : 
00741 //=======================================================================
00743 void GRID::getCellPosition(const int Number, int& i, int& j, int& k) const
00744   throw (MEDEXCEPTION)
00745 {
00746   
00747   const char* LOC = "GRID::getCellPosition(Number, i,j,k) :";
00748   BEGIN_OF_MED(LOC);
00749 
00750   int Len[4] = {0,_iArrayLength-1, _jArrayLength-1, _kArrayLength-1 }, I=1, J=2;
00751   // , K=3; !! UNUSED VARIABLE !!
00752 
00753 //  if (Number <= 0 || Number > getCellNumber(Len[I]-1, Len[J]-1, Len[K]-1))
00754 //    throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
00755 
00756   int ijLen = Len[I] * Len[J]; // nb in a full k layer
00757   int kLen = (Number - 1) % ijLen;    // nb in the non full k layer
00758   
00759   i = kLen % Len[J];
00760   j = kLen / Len[J];
00761   k = (Number - 1) / ijLen;
00762 
00763   END_OF_MED(LOC);
00764 }
00765 
00766 //=======================================================================
00767 //function : getEdgePosition
00768 //purpose  : 
00769 //=======================================================================
00771 void GRID::getEdgePosition(const int Number, int& Axis, int& i, int& j, int& k)
00772   const throw (MEDEXCEPTION)
00773 {
00774   const char * LOC = "GRID::getEdgePosition(Number, i,j,k) :";
00775 
00776   BEGIN_OF_MED(LOC);
00777 
00778   if (!_jArrayLength)
00779     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "no edges in the grid: "));
00780   
00781   if (Number <= 0)
00782     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
00783 
00784   int Len[4] = {0,_iArrayLength, _jArrayLength, _kArrayLength }, I=1, J=2, K=3;
00785 
00786   int theNb = Number;
00787   int maxAxis = _kArrayLength ? 3 : 2;
00788   
00789   for (Axis = 1; Axis <= maxAxis; ++Axis)
00790   {
00791     Len[Axis]--;  // one less edge in Axis direction
00792 
00793     // max edge number in Axis direction
00794     int maxNb = getEdgeNumber (Axis, Len[I]-1, Len[J]-1, Len[K]-1);
00795     
00796     if (theNb > maxNb)
00797     {
00798       Len[Axis]++;
00799       theNb -= maxNb;
00800       continue;
00801     }
00802     
00803     if (theNb == maxNb)
00804     {
00805       i = Len[I]-1;
00806       j = Len[J]-1;
00807       k = Len[K]-1;
00808     }
00809     else
00810     {
00811       int ijLen = Len[I] * Len[J]; // nb in a full k layer
00812       int kLen = (theNb - 1) % ijLen;    // nb in the non full k layer
00813       i = kLen % Len[J];
00814       j = kLen / Len[J];
00815       k = (theNb - 1) / ijLen;
00816     }
00817 
00818   END_OF_MED(LOC);
00819 
00820     return;
00821   }
00822   
00823   throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
00824 }
00825 
00826 //=======================================================================
00827 //function : getFacePosition
00828 //purpose  : return position (i,j,k) of an entity #Number
00829 //           Axis [1,2,3] means one of directions: along i, j or k
00830 //           For Cell contituents (FACE or EDGE), Axis selects one of those having same (i,j,k):
00831 //           * a FACE which is normal to direction along given Axis;
00832 //           * an EDGE going along given Axis.
00833 //           Exception for Number out of range
00834 //=======================================================================
00836 void GRID::getFacePosition(const int Number, int& Axis, int& i, int& j, int& k)
00837   const throw (MEDEXCEPTION)
00838 {
00839   const char * LOC = "GRID::getFacePosition(Number, i,j,k) :";
00840 
00841   BEGIN_OF_MED(LOC);
00842 
00843   if (_kArrayLength == 0) {
00844     getCellPosition(Number, i, j, k);
00845     Axis = 1;
00846     return;
00847   };
00848 
00849   if (!_kArrayLength)
00850     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "no faces in the grid: "));
00851   
00852   if ( Number <= 0 )
00853     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
00854 
00855   int Len[4] = {0,_iArrayLength-1, _jArrayLength-1, _kArrayLength-1 }, I=1, J=2, K=3;
00856   int theNb = Number;
00857   
00858   for (Axis = 1; Axis <= 3; ++Axis)
00859   {
00860     Len[Axis]++;
00861 
00862     // max face number in Axis direction
00863     int maxNb = getFaceNumber (Axis, Len[I]-1, Len[J]-1, Len[K]-1);
00864     
00865     if (theNb > maxNb)
00866     {
00867       Len[Axis]--;
00868       theNb -= maxNb;
00869       continue;
00870     }
00871     
00872     if (theNb == maxNb)
00873     {
00874       i = Len[I]-1;
00875       j = Len[J]-1;
00876       k = Len[K]-1;
00877     }
00878     else
00879     {
00880       int ijLen = Len[I] * Len[J]; // nb in a full k layer
00881       int kLen = (theNb - 1)  % ijLen;    // nb in the non full k layer
00882       i = kLen % Len[J];
00883       j = kLen / Len[J];
00884       k = (theNb - 1) / ijLen;
00885     }
00886 
00887   END_OF_MED(LOC);
00888 
00889     return;
00890   }
00891   
00892   throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
00893 }
00901 //================================================================================
00908 //================================================================================
00909 
00910 int GRID::getNumberOfTypes(MED_EN::medEntityMesh entity) const
00911 {
00912   MESSAGE_MED("GRID::getNumberOfTypes(medEntityMesh entity) : "<<entity);
00913   return 1; // a grid has one type
00914 }
00915 
00916 //================================================================================
00920 //================================================================================
00921 
00922 int GRID::getNumberOfElements(MED_EN::medEntityMesh entity, MED_EN::medGeometryElement Type) const
00923 {
00924   int numberOfElements=0;
00925 
00926   // Cas où le nombre d'éléments n'est pas nul
00927   if (entity==MED_EN::MED_FACE && (Type==MED_EN::MED_QUAD4 || Type==MED_EN::MED_ALL_ELEMENTS) && getMeshDimension()>2)
00928     numberOfElements=
00929       (_iArrayLength-1)*(_jArrayLength-1)*(_kArrayLength  )+
00930       (_iArrayLength-1)*(_jArrayLength  )*(_kArrayLength-1)+
00931       (_iArrayLength  )*(_jArrayLength-1)*(_kArrayLength-1);
00932 
00933   else if (entity==MED_EN::MED_EDGE && (Type==MED_EN::MED_SEG2 || Type==MED_EN::MED_ALL_ELEMENTS))
00934     if ( _spaceDimension==2)
00935       numberOfElements=_iArrayLength*(_jArrayLength-1) + (_iArrayLength-1)*_jArrayLength;
00936     else if ( _spaceDimension==1)
00937       numberOfElements=_iArrayLength-1;
00938     else // 3D
00939       numberOfElements=
00940         (_iArrayLength*(_jArrayLength-1) + (_iArrayLength-1)*_jArrayLength) * _kArrayLength +
00941         _iArrayLength*_jArrayLength*(_kArrayLength-1);
00942 
00943   else if (entity==MED_EN::MED_NODE && (Type==MED_EN::MED_NONE || Type==MED_EN::MED_ALL_ELEMENTS) && _spaceDimension>0)
00944     numberOfElements=getNumberOfNodes();
00945 
00946   else if (entity==MED_EN::MED_CELL && _spaceDimension==3 && (Type==MED_EN::MED_HEXA8 || Type==MED_EN::MED_ALL_ELEMENTS) )
00947     numberOfElements=(_iArrayLength-1)*(_jArrayLength-1)*(_kArrayLength-1);
00948 
00949   else if (entity==MED_EN::MED_CELL && _spaceDimension==2 && (Type==MED_EN::MED_QUAD4 || Type==MED_EN::MED_ALL_ELEMENTS))
00950     numberOfElements=(_iArrayLength-1)*(_jArrayLength-1);
00951 
00952   else if (entity==MED_EN::MED_CELL && _spaceDimension==1 && (Type==MED_EN::MED_SEG2 || Type==MED_EN::MED_ALL_ELEMENTS) )
00953     numberOfElements=_iArrayLength-1;
00954 
00955   MESSAGE_MED("GRID::getNumberOfElements - entity=" << entity << " Type=" << Type);
00956   MESSAGE_MED("_spaceDimension=" << _spaceDimension << "  numberOfElements=" << numberOfElements);
00957 
00958   return numberOfElements;
00959 }
00960 
00961 //================================================================================
00965 //================================================================================
00966 
00967 MED_EN::medGeometryElement GRID::getElementType(MED_EN::medEntityMesh Entity,int Number) const
00968 {
00969   return getTypes(Entity)[0];
00970 }
00971 
00972 //================================================================================
00976 //================================================================================
00977 
00978 int GRID::getMeshDimension() const
00979 {
00980   return getSpaceDimension();
00981 }
00982 
00983 //================================================================================
00987 //================================================================================
00988 
00989 bool GRID::getIsAGrid() const
00990 {
00991   return true;
00992 }
00993 
00994 //================================================================================
00998 //================================================================================
00999 
01000 int GRID::getNumberOfNodes() const
01001 {
01002   if ( _gridType == MED_EN::MED_BODY_FITTED )
01003     return _coordinate ? _coordinate->getNumberOfNodes() : 0;
01004 
01005   switch ( _spaceDimension )
01006   {
01007   case 3: return _iArrayLength * _jArrayLength * _kArrayLength;
01008   case 2: return _iArrayLength * _jArrayLength;
01009   case 1: return _iArrayLength;
01010   }
01011   return 0;
01012 }
01013 
01014 //=======================================================================
01018 //=======================================================================
01019 
01020 std::string         GRID::getCoordinatesSystem() const
01021 {
01022   return _coordinate ? _coordinate->getCoordinatesSystem() : defaultStrings()[0];
01023 }
01024 
01025 //=======================================================================
01033 //=======================================================================
01034 
01035 const std::string * GRID::getCoordinatesNames() const
01036 {
01037   return _coordinate ? _coordinate->getCoordinatesNames() : defaultStrings();
01038 }
01039 
01040 //=======================================================================
01045 //=======================================================================
01046 
01047 const std::string * GRID::getCoordinatesUnits() const
01048 {
01049   return _coordinate ? _coordinate->getCoordinatesUnits() : defaultStrings();
01050 }
01051 
01052 //=======================================================================
01066 //=======================================================================
01067 
01068 SUPPORT * GRID::getBoundaryElements(MED_EN::medEntityMesh Entity) const throw (MEDEXCEPTION)
01069 {
01070   const char * LOC = "GMESH::getBoundaryElements() : " ;
01071   if ( _spaceDimension < 2 )
01072     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not implemented in " << _spaceDimension <<"D space !"));
01073 
01074   if ( _gridType == MED_POLAR) 
01075     throw MEDEXCEPTION("GRID::getBoundaryElements() : not implemented on MED_POLAR grig");
01076 
01077 
01078   medEntityMesh entityToParse=Entity;
01079   if(_spaceDimension == 3)
01080     if (Entity != MED_FACE)
01081       {
01082         if(Entity==MED_NODE)
01083           entityToParse=MED_FACE;
01084         else
01085           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<Entity<<" !"));
01086       }
01087   if(_spaceDimension == 2)
01088     if(Entity != MED_EDGE)
01089       {
01090         if(Entity==MED_NODE)
01091           entityToParse=MED_EDGE;
01092         else
01093           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<Entity<<" !"));
01094       }
01095   list<int> bnd_elems1, bnd_elems2;
01096   int numberOf = getNumberOfElements(entityToParse,MED_ALL_ELEMENTS) ;
01097 
01098   if ( _coordinate->getNumberOfNodes() > 0 ) // BODY FITTED
01099     {
01100       throw MEDEXCEPTION("GRID::getBoundaryElements() : not implemented on BOBY FITTED grig");
01101     }
01102   else if ( entityToParse == MED_FACE ) // 3D CARTESIAN
01103     {
01104       const int nb_x = getArrayLength( 1 ) - 1;
01105       const int nb_y = getArrayLength( 2 ) - 1;
01106       const int nb_z = getArrayLength( 3 ) - 1;
01107       // normal to OX
01108       for ( int z = 0; z < nb_z; ++z )
01109         for ( int y = 0; y < nb_y; ++y )
01110           {
01111             bnd_elems1.push_back( getFaceNumber( 1, 0,    y, z ));
01112             bnd_elems2.push_back( getFaceNumber( 1, nb_x, y, z ));
01113           }
01114       bnd_elems1.splice( bnd_elems1.end(), bnd_elems2 ); // to have ids in increasing order
01115 
01116       // normal to OY
01117       for ( int z = 0; z < nb_z; ++z )
01118         for ( int x = 0; x < nb_x; ++x )
01119           {
01120             bnd_elems1.push_back( getFaceNumber( 2, x, 0,    z ));
01121             bnd_elems2.push_back( getFaceNumber( 2, x, nb_y, z ));
01122           }
01123       bnd_elems1.splice( bnd_elems1.end(), bnd_elems2 );
01124 
01125       // normal to OZ
01126       for ( int y = 0; y < nb_y; ++y )
01127         for ( int x = 0; x < nb_x; ++x )
01128           {
01129             bnd_elems1.push_back( getFaceNumber( 3, x, y, 0    ));
01130             bnd_elems2.push_back( getFaceNumber( 3, x, y, nb_z ));
01131           }
01132       bnd_elems1.splice( bnd_elems1.end(), bnd_elems2 );
01133     }
01134   else
01135     {
01136       const int nb_x = getArrayLength( 1 ) - 1;
01137       const int nb_y = getArrayLength( 2 ) - 1;
01138       // edge || OX
01139       for ( int x = 0; x < nb_x; ++x )
01140         {
01141           bnd_elems1.push_back( getEdgeNumber( 1, x, 0 ));
01142           bnd_elems2.push_back( getEdgeNumber( 1, x, nb_y ));
01143         }
01144       bnd_elems1.splice( bnd_elems1.end(), bnd_elems2 ); // to have ids in increasing order
01145       // edge || OY
01146       for ( int y = 0; y < nb_y; ++y )
01147         {
01148           bnd_elems1.push_back( getEdgeNumber( 2, y, 0 ));
01149           bnd_elems2.push_back( getEdgeNumber( 2, y, nb_x ));
01150         }
01151       bnd_elems1.splice( bnd_elems1.end(), bnd_elems2 );
01152     }
01153 
01154   if ( bnd_elems1.empty() && numberOf != 0 )
01155     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No boundary elements found by reverse descending connectivity for entity "<<Entity<<" !"));
01156 
01157   if ( Entity == MED_NODE )
01158     return buildSupportOnNodeFromElementList(bnd_elems1,entityToParse);
01159   else
01160     return buildSupportOnElementsFromElementList(bnd_elems1,entityToParse);
01161 }
01162 
01163 SUPPORT * GRID::getSkin(const SUPPORT * Support3D) throw (MEDEXCEPTION)
01164 {
01165   throw MEDEXCEPTION("GRID::getSkin() : Not implemented yet");
01166 }
01167 SUPPORT *GRID::buildSupportOnNodeFromElementList(const std::list<int>& listOfElt,
01168                                                  MED_EN::medEntityMesh entity) const
01169   throw (MEDEXCEPTION)
01170 {
01171   throw MEDEXCEPTION("GRID::buildSupportOnNodeFromElementList() : Not implemented yet");
01172 }
01173 void GRID::fillSupportOnNodeFromElementList(const std::list<int>& listOfElt,
01174                                             SUPPORT *             supportToFill) const
01175   throw (MEDEXCEPTION)
01176 {
01177   throw MEDEXCEPTION("GRID::fillSupportOnNodeFromElementList() : Not implemented yet");
01178 }
01179 
01180 FIELD<double>* GRID::getVolume (const SUPPORT * Support, bool isAbs) const
01181   throw (MEDEXCEPTION)
01182 {
01183   throw MEDEXCEPTION("GRID::getVolume() : Not implemented yet");
01184 }
01185 
01186 FIELD<double>* GRID::getArea (const SUPPORT * Support) const
01187   throw (MEDEXCEPTION)
01188 {
01189   throw MEDEXCEPTION("GRID::getArea() : Not implemented yet");
01190 }
01191 
01192 FIELD<double>* GRID::getLength (const SUPPORT * Support) const
01193   throw (MEDEXCEPTION)
01194 {
01195   throw MEDEXCEPTION("GRID::getLength() : Not implemented yet");
01196 }
01197 
01198 FIELD<double>* GRID::getNormal (const SUPPORT * Support) const
01199   throw (MEDEXCEPTION)
01200 {
01201   throw MEDEXCEPTION("GRID::getNormal() : Not implemented yet");
01202 }
01203 
01204 FIELD<double>* GRID::getBarycenter (const SUPPORT * Support) const
01205   throw (MEDEXCEPTION)
01206 {
01207   throw MEDEXCEPTION("GRID::getBarycenter() : Not implemented yet");
01208 }
01209 
01210 vector< vector<double> >   GRID::getBoundingBox() const
01211 {
01212   throw MEDEXCEPTION("GRID::getBoundingBox() : Not implemented yet");
01213 }