src/MEDMEM/MEDMEM_Grid.cxx

   1
   2 // Copyright (C) 2005  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   3 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   4 //
   5 // This library is free software; you can redistribute it and/or
   6 // modify it under the terms of the GNU Lesser General Public
   7 // License as published by the Free Software Foundation; either
   8 // version 2.1 of the License.
   9 //
  10 // This library is distributed in the hope that it will be useful
  11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 // Lesser General Public License for more details.
  14 //
  15 // You should have received a copy of the GNU Lesser General Public
  16 // License along with this library; if not, write to the Free Software
  17 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 //
  19 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 //
  21 // File      : MEDMEM_Grid.hxx
  22 // Created   : Wed Dec 18 08:35:26 2002
  23 // Descr     : class containing structured mesh data
  24
  25 // Author    : Edward AGAPOV (eap)
  26 // Project   : SALOME Pro
  27 // Module    : MED
  28 // $Header$
  29
  30
  31 #include "MEDMEM_Grid.hxx"
  32 #include "MEDMEM_CellModel.hxx"
  33 #include "MEDMEM_SkyLineArray.hxx"
  34 #include "MEDMEM_DriverFactory.hxx"
  35
  36 using namespace std;
  37 using namespace MEDMEM;
  38 using namespace MED_EN;
  39
  40 //=======================================================================
  41 //function : GRID
  42 //purpose  : empty constructor
  43 //=======================================================================
  44
  45 GRID::GRID() {
  46   init();
  47   MESSAGE("A GRID CREATED");
  48 }
  49 //
  50 //=======================================================================
  51 ////function : GRID
  52 ////purpose  : array constructor
  53 ////=======================================================================
  54 GRID::GRID(const std::vector<std::vector<double> >& xyz_array,
  55            const std::vector<std::string>&          coord_name,
  56            const std::vector<std::string>&          coord_unit,
  57            const med_grid_type                      type) : _gridType(type)
  58 {
  59     init(); // PAL 12136
  60     _is_default_gridType = false;
  61
  62     _spaceDimension = xyz_array.size();
  63     _meshDimension  = _spaceDimension; // PAL14482
  64
  65     // compute & set _numberOfNodes
  66     int NumberOfNodes=1 ;
  67     for(int i=0;i!=xyz_array.size();++i)
  68         NumberOfNodes*=xyz_array[i].size();
  69     _numberOfNodes = NumberOfNodes ;
  70
  71     // create a non allocated COORDINATE
  72     _coordinate = new COORDINATE(_spaceDimension, &coord_name[0], &coord_unit[0]);
  73     string coordinateSystem = "UNDEFINED";
  74     if( _gridType == MED_CARTESIAN)
  75         coordinateSystem = "CARTESIAN";
  76     else if ( _gridType == MED_POLAR)
  77         coordinateSystem = "CYLINDRICAL";
  78     _coordinate->setCoordinatesSystem(coordinateSystem);
  79
  80     // set the GRID part
  81     if (_spaceDimension>=1)
  82     {
  83         _iArrayLength=xyz_array[0].size();
  84         _iArray=new double[_iArrayLength];
  85         std::copy(xyz_array[0].begin(),xyz_array[0].end(),_iArray);
  86     }
  87     if (_spaceDimension>=2)
  88     {
  89         _jArrayLength=xyz_array[1].size();
  90         _jArray=new double[_jArrayLength];
  91         std::copy(xyz_array[1].begin(),xyz_array[1].end(),_jArray);
  92     }
  93     if (_spaceDimension>=3)
  94     {
  95         _kArrayLength=xyz_array[2].size();
  96         _kArray=new double[_kArrayLength];
  97         std::copy(xyz_array[2].begin(),xyz_array[2].end(),_kArray);
  98     }
  99
 100     _is_coordinates_filled  = false;
 101     _is_connectivity_filled = false;
 102     _isAGrid = true;
 103 }
 104
 105 //=======================================================================
 106 //function : GRID
 107 //purpose  : empty constructor
 108 //=======================================================================
 109
 110 GRID::GRID(const med_grid_type type)
 111 {
 112   init();
 113   _gridType = type;
 114   _is_default_gridType = false;
 115   MESSAGE("A TYPED GRID CREATED");
 116 }
 117
 118 //=======================================================================
 119 //function : GRID
 120 //purpose  : copying constructor
 121 //=======================================================================
 122
 123 GRID::GRID(const GRID& otherGrid) {
 124   *this = otherGrid;
 125 }
 126
 127 //=======================================================================
 128 //function : ~GRID
 129 //purpose  : destructor
 130 //=======================================================================
 131
 132 GRID::~GRID() {
 133   MESSAGE("GRID::~GRID() : Destroying the Grid");
 134   if ( _iArray != (double* ) NULL) delete [] _iArray;
 135   if ( _jArray != (double* ) NULL) delete [] _jArray;
 136   if ( _kArray != (double* ) NULL) delete [] _kArray;
 137 }
 138
 139 //=======================================================================
 140 //function : init
 141 //purpose :
 142 //=======================================================================
 143
 144 void GRID::init()
 145 {
 146   _gridType = MED_CARTESIAN;
 147   _is_default_gridType = true;
 148
 149   _iArray = _jArray = _kArray = (double* ) NULL;
 150   _iArrayLength = _jArrayLength = _kArrayLength = 0;
 151
 152   _is_coordinates_filled  = false;
 153   _is_connectivity_filled = false;
 154
 155   MESH::init();
 156
 157   _isAGrid = true;
 158 }
 159
 160 //=======================================================================
 161 //function: operator=
 162 //purpose : operator=
 163 //=======================================================================
 164
 165 GRID & GRID::operator=(const GRID & otherGrid)
 166 {
 167   // NOT IMPLEMENTED
 168
 169   MESH::operator=(otherGrid);
 170   return *this;
 171 }
 172
 173 //=======================================================================
 174 //function : GRID
 175 //purpose  : Create a GRID object using a MESH driver of type
 176 //           (MED_DRIVER, ....) associated with file <fileName>.
 177 //=======================================================================
 178
 179 GRID::GRID(driverTypes driverType, const string &  fileName,
 180            const string &  driverName)// : MESH(driverType, fileName, driverName)
 181 {
 182   const char * LOC ="GRID::GRID(driverTypes , const string &  , const string &) : ";
 183
 184   BEGIN_OF(LOC);
 185
 186   init();
 187   GENDRIVER *myDriver=DRIVERFACTORY::buildDriverForMesh(driverType,fileName,this,driverName,MED_LECT);
 188   int current = addDriver(*myDriver);
 189   delete myDriver;
 190   _drivers[current]->open();
 191   _drivers[current]->read();
 192   _drivers[current]->close();
 193
 194   fillMeshAfterRead();
 195
 196   END_OF(LOC);
 197 };
 198
 199 /*!
 200   return the GRID Geometric type, without computing all connectivity
 201 */
 202 const medGeometryElement * GRID::getTypes(medEntityMesh entity) const
 203 {
 204     static const medGeometryElement _gridGeometry[4]={MED_HEXA8,MED_QUAD4,MED_SEG2,MED_POINT1};
 205     int i=0;
 206     if(entity==MED_CELL)
 207     {
 208         i=3-_spaceDimension;
 209     }
 210     else if(entity==MED_FACE && _spaceDimension>2 )
 211         i=1;
 212     else if(entity==MED_EDGE && _spaceDimension>1 )
 213         i=2;
 214     else if(entity==MED_NODE && _spaceDimension>0)
 215         i=3;
 216     else
 217         throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::getGeometricTypes : Entity not defined !"));
 218     return &_gridGeometry[i];
 219 }
 220
 221 MED_EN::medGeometryElement * GRID::getTypesWithPoly(MED_EN::medEntityMesh Entity) const
 222 {
 223   int size=getNumberOfTypesWithPoly(Entity);
 224   MED_EN::medGeometryElement *ret=new MED_EN::medGeometryElement[size];
 225   memcpy(ret,getTypes(Entity),size*sizeof(MED_EN::medGeometryElement));
 226   return ret;
 227 }
 228
 229 //=======================================================================
 230 //function : fillMeshAfterRead
 231 //purpose  :
 232 //=======================================================================
 233
 234 void GRID::fillMeshAfterRead()
 235 {
 236  // fill not only MESH  (:-)
 237   _is_coordinates_filled  = false;
 238   _is_connectivity_filled = false;
 239
 240   if (_gridType == MED_BODY_FITTED)
 241   {
 242     _is_coordinates_filled = true;
 243
 244     // nb of nodes in each direction is not known, set anything
 245     // in order to be able to work anyhow
 246 //     INFOS("GRID::fillMeshAfterRead(): This stub must be removed");
 247 //     switch (_spaceDimension) {
 248 //     case 1:
 249 //       _iArrayLength = _numberOfNodes;
 250 //       _jArrayLength = 0;
 251 //       _kArrayLength = 0;
 252 //       break;
 253 //     case 2:
 254 //       _iArrayLength = 2;
 255 //       _jArrayLength = _numberOfNodes / _iArrayLength;
 256 //       _kArrayLength = 0;
 257 //       break;
 258 //     default:
 259 //       _iArrayLength = 2;
 260 //       _jArrayLength = _numberOfNodes / _iArrayLength / 2;
 261 //       _kArrayLength = _numberOfNodes / _iArrayLength / _jArrayLength;
 262 //     }
 263     //cout << "ARRAY LENGTHS: " << _iArrayLength << " " << _jArrayLength
 264 //      << " " << _kArrayLength << endl;
 265   }
 266 //   else {
 267 //     int NbNodes = _iArrayLength;
 268 //     if (_jArrayLength)
 269 //       NbNodes *= _jArrayLength;
 270 //     if (_kArrayLength)
 271 //       NbNodes *= _kArrayLength;
 272 //     MESH::_numberOfNodes = NbNodes;
 273 //   }
 274
 275   MESH::_meshDimension = MESH::_spaceDimension;
 276 }
 277
 278 //=======================================================================
 279 //function : fillCoordinates
 280 //purpose  :
 281 //=======================================================================
 282
 283 void GRID::fillCoordinates() const
 284 {
 285   if (_is_coordinates_filled)
 286     return;
 287
 288   const char * LOC ="GRID::fillCoordinates()";
 289   BEGIN_OF(LOC);
 290
 291   // if coordonate has not been allocated, perform shalow copy, transfer ownership of matrix
 292   if(_coordinate->getSpaceDimension()*_coordinate->getNumberOfNodes() == 0)
 293       _coordinate->setCoordinates(new MEDARRAY<double>(_spaceDimension,_numberOfNodes,MED_FULL_INTERLACE),true);
 294
 295   double* myCoord = const_cast <double *> ( _coordinate->getCoordinates(MED_FULL_INTERLACE) );
 296
 297   bool hasJ = _jArrayLength, hasK = _kArrayLength;
 298   int J = hasJ ? _jArrayLength : 1;
 299   int K = hasK ? _kArrayLength : 1;
 300   //int nb, !! UNUSED VARIABLE !!
 301   int i, j, k;
 302   for (k=0; k < K; ++k) {
 303     for (j=0; j < J; ++j) {
 304       for (i=0; i < _iArrayLength; ++i) {
 305
 306         * myCoord = _iArray[ i ];
 307         ++ myCoord;
 308
 309         if (hasJ)
 310         {
 311           * myCoord = _jArray[ j ];
 312           ++ myCoord;
 313
 314           if (hasK)
 315           {
 316             * myCoord = _kArray[ k ];
 317             ++ myCoord;
 318           }
 319         }
 320       }
 321     }
 322   }
 323
 324   (const_cast <GRID *> (this))->_is_coordinates_filled = true;
 325   END_OF(LOC);
 326 }
 327
 328 //=======================================================================
 329 //function : makeConnectivity
 330 //purpose  :
 331 //=======================================================================
 332
 333 CONNECTIVITY * GRID::makeConnectivity (medEntityMesh           Entity,
 334                                        medGeometryElement Geometry,
 335                                        int                NbEntities,
 336                                        int                NbNodes,
 337                                        int                nbMeshNodes,
 338                                        const int *                    NodeNumbers)
 339   const
 340 {
 341   CONNECTIVITY * Connectivity     = new CONNECTIVITY(Entity) ;
 342   Connectivity->_numberOfNodes    = nbMeshNodes ;
 343   Connectivity->_entityDimension  = Geometry/100 ;
 344
 345   int numberOfGeometricType    = 1;
 346   Connectivity->_numberOfTypes = numberOfGeometricType;
 347
 348   Connectivity->_count    = new int [numberOfGeometricType + 1] ;
 349   Connectivity->_count[0] = 1;
 350   Connectivity->_count[1] = 1 + NbEntities;
 351
 352   Connectivity->_type    = new CELLMODEL [numberOfGeometricType];
 353   Connectivity->_type[0] = CELLMODEL( Geometry ) ;
 354
 355   Connectivity->_geometricTypes    = new medGeometryElement [ numberOfGeometricType ];
 356   Connectivity->_geometricTypes[0] = Geometry;
 357
 358   //  Connectivity->_nodal = new MEDSKYLINEARRAY() ;
 359   int * skyLineArrayIndex = new int [NbEntities + 1];
 360   int i, j, nbEntityNodes = Connectivity->_type[0].getNumberOfNodes();
 361   for (i=0, j=1; i <= NbEntities; ++i, j += nbEntityNodes)
 362     skyLineArrayIndex [ i ] = j;
 363
 364   //  Connectivity->_nodal->setMEDSKYLINEARRAY (NbEntities, NbNodes,
 365   //skyLineArrayIndex, NodeNumbers);
 366
 367   Connectivity->_nodal = new MEDSKYLINEARRAY (NbEntities, NbNodes,
 368                                                skyLineArrayIndex, NodeNumbers);
 369
 370   delete [] skyLineArrayIndex;
 371
 372   // test connectivity right away
 373 //   for (med_int j = 0; j < numberOfGeometricType; j++)
 374 //   {
 375 //     int node_number = Connectivity->_type[j].getNumberOfNodes();
 376 //     for (med_int k = Connectivity->_count[j]; k < Connectivity->_count[j+1]; k++)
 377 //       for (med_int local_node_number = 1 ; local_node_number < node_number+1; local_node_number++)
 378 //       {
 379 //         cout << "MEDSKYLINEARRAY::getIJ(" << k << ", " << local_node_number << endl;
 380 //         med_int global_node = Connectivity->_nodal->getIJ(k,local_node_number) ;
 381 //         cout << "...= " << global_node << endl;
 382 //       }
 383 //   }
 384
 385   return Connectivity;
 386 }
 387
 388 //=======================================================================
 389 //function : fillConnectivity
 390 //purpose  : fill _coordinates and _connectivity of MESH if not yet done
 391 //=======================================================================
 392
 393 void GRID::fillConnectivity() const
 394 {
 395   if (_is_connectivity_filled)
 396   {
 397     MESSAGE("GRID::fillConnectivity(): Already filled");
 398     return;
 399   }
 400
 401   const char * LOC = "GRID::fillConnectivity() ";
 402   BEGIN_OF(LOC);
 403
 404   int nbCells, nbFaces, nbEdges;
 405   int nbCNodes, nbFNodes, nbENodes, nbMeshNodes;
 406   int indexC, indexF, indexE;
 407   int * nodeCNumbers, * nodeFNumbers, * nodeENumbers;
 408   // about descending connectivity
 409   int nbSub, nbRevSub, indexSub, indexRevSub, * subNumbers, * subRevNumbers;
 410
 411   bool hasFaces = _kArrayLength, hasEdges = _jArrayLength;
 412
 413   int iLenMin1 = _iArrayLength-1, jLenMin1 = _jArrayLength-1;
 414   int kLenMin1 = _kArrayLength-1, ijLenMin1 = iLenMin1 * jLenMin1;
 415   if (!hasEdges) jLenMin1 = 1;
 416   if (!hasFaces) kLenMin1 = 1;
 417
 418   // nb of cells and of their connectivity nodes
 419
 420   nbCells = iLenMin1 * jLenMin1 * kLenMin1;
 421   nbMeshNodes = _iArrayLength * (_jArrayLength ? _jArrayLength : 1) * (_kArrayLength ? _kArrayLength : 1);
 422   nbCNodes = nbCells * 2 * (hasEdges ? 2 : 1) * (hasFaces ? 2 : 1);
 423   nodeCNumbers = new int [ nbCNodes ];
 424
 425   // nb of faces and of their connectivity nodes
 426
 427   if (hasFaces) {
 428     nbFaces  = _iArrayLength * jLenMin1 * kLenMin1;
 429     nbFaces += _jArrayLength * kLenMin1 * iLenMin1;
 430     nbFaces += _kArrayLength * iLenMin1 * jLenMin1;
 431     nbFNodes = nbFaces * 4;
 432     nodeFNumbers = new int [ nbFNodes ];
 433   } else
 434     nbFaces = nbFNodes = 0;
 435
 436   // nb of edges and of their connectivity nodes
 437
 438   if (hasEdges) {
 439     if (_kArrayLength) { // 3d grid
 440       nbEdges  = iLenMin1 * _jArrayLength * _kArrayLength;
 441       nbEdges += jLenMin1 * _kArrayLength * _iArrayLength;
 442       nbEdges += kLenMin1 * _iArrayLength * _jArrayLength;
 443     }
 444     else if (_jArrayLength) { // 2d
 445       nbEdges  = iLenMin1 * _jArrayLength;
 446       nbEdges += jLenMin1 * _iArrayLength;
 447     }
 448     nbENodes = nbEdges * 2;
 449     nodeENumbers = new int [ nbENodes ];
 450   } else
 451     nbEdges = nbENodes = 0;
 452
 453   // nb of descending connectivity Elements
 454
 455   if (hasFaces)
 456   {
 457     nbSub = nbCells * 6;
 458     nbRevSub = nbFaces * 2;
 459   }
 460   else
 461   {
 462     nbSub = nbCells * 4;
 463     nbRevSub = nbEdges * 2;
 464   }
 465   subNumbers = new int [ nbSub ];
 466   subRevNumbers = new int [ nbRevSub ];
 467   for (int r=0; r<nbRevSub; ++r)
 468     subRevNumbers[ r ] = 0;
 469
 470   int nSubI = 1, nSubJ, nSubK; // subelement numbers
 471   if (hasFaces)
 472   {
 473     nSubJ = getFaceNumber(2, 0, 0, 0);
 474     nSubK = getFaceNumber(3, 0, 0, 0);
 475   }
 476   else
 477     nSubJ = getEdgeNumber(2, 0, 0, 0);
 478
 479
 480   // fill cell node numbers and descending element numbers
 481
 482
 483   indexC = indexF = indexE = indexSub = indexRevSub = -1;
 484   int iNode = 0, iCell = 0;
 485   int ijLen = _iArrayLength * _jArrayLength;
 486   int i, j, k, n1, n2, n3 ,n4;
 487
 488   int I = _iArrayLength;
 489   int J = _jArrayLength ? _jArrayLength : 2; // pass at least once
 490   int K = _kArrayLength ? _kArrayLength : 2;
 491   // pass by all but last nodes in all directions
 492   for (k = 1; k < K; ++k ) {
 493     for (j = 1; j < J; ++j ) {
 494       for (i = 1; i < I; ++i ) {
 495
 496         ++iCell;
 497
 498         n1 = ++iNode; // iNode
 499         n2 = n1 + 1;
 500         nodeCNumbers [ ++indexC ] = n1;
 501         nodeCNumbers [ ++indexC ] = n2;
 502
 503         if (hasEdges) { // at least 2d
 504           n3 = n2 + I;
 505           n4 = n3 - 1;
 506           nodeCNumbers [ ++indexC ] = n3;
 507           nodeCNumbers [ ++indexC ] = n4;
 508         }
 509         if (hasFaces) { // 3d
 510           nodeCNumbers [ ++indexC ] = n1 + ijLen;
 511           nodeCNumbers [ ++indexC ] = n2 + ijLen;
 512           nodeCNumbers [ ++indexC ] = n3 + ijLen;
 513           nodeCNumbers [ ++indexC ] = n4 + ijLen;
 514
 515           // descending faces
 516           n1 = nSubI;
 517           n2 = n1 + 1;
 518           n3 = (n1-1) * 2;
 519           n4 = (n2-1) * 2 + 1;
 520           subNumbers [ ++indexSub ] = n1;
 521           subRevNumbers[ n3 ] = iCell;
 522           subNumbers [ ++indexSub ] = n2;
 523           subRevNumbers[ n4 ] = -iCell;
 524           n1 = nSubJ;
 525           n2 = n1 + iLenMin1;
 526           n3 = (n1-1) * 2;
 527           n4 = (n2-1) * 2 + 1;
 528           subNumbers [ ++indexSub ] = n1;
 529           subRevNumbers[ n3 ] = iCell;
 530           subNumbers [ ++indexSub ] = n2;
 531           subRevNumbers[ n4 ] = -iCell;
 532           n1 = nSubK;
 533           n2 = n1 + ijLenMin1;
 534           n3 = (n1-1) * 2;
 535           n4 = (n2-1) * 2 + 1;
 536           subNumbers [ ++indexSub ] = n1;
 537           subRevNumbers[ n3 ] = iCell;
 538           subNumbers [ ++indexSub ] = n2;
 539           subRevNumbers[ n4 ] = -iCell;
 540         }
 541         else
 542         {
 543           // descending edges
 544           n1 = nSubI;
 545           n2 = n1 + iLenMin1;
 546           n3 = (n1-1) * 2;
 547           n4 = (n2-1) * 2 + 1;
 548           subNumbers [ ++indexSub ] = n1;
 549           subRevNumbers[ n3 ] = iCell;
 550           subNumbers [ ++indexSub ] = n2;
 551           subRevNumbers[ n4 ] = -iCell;
 552           n1 = nSubJ;
 553           n2 = n1 + 1;
 554           n3 = (n1-1) * 2;
 555           n4 = (n2-1) * 2 + 1;
 556           subNumbers [ ++indexSub ] = n1;
 557           subRevNumbers[ n3 ] = iCell;
 558           subNumbers [ ++indexSub ] = n2;
 559           subRevNumbers[ n4 ] = -iCell;
 560         }
 561         ++nSubI; ++nSubJ; ++nSubK;
 562       }
 563       ++iNode; // skip the last node in a row
 564       if (hasFaces)
 565         ++nSubI;
 566       else
 567         ++nSubJ;
 568     }
 569     iNode += I; // skip the whole last row
 570   }
 571
 572   // fill face node numbers
 573
 574   int ax, AX = hasFaces ? 3 : 0;
 575   for (ax = 1; ax <= AX; ++ax) {
 576
 577     iNode = 0;
 578
 579     I = _iArrayLength;
 580     J = _jArrayLength;
 581     K = _kArrayLength;
 582     switch (ax) {
 583     case 1:  --J; --K; break;
 584     case 2:  --K; --I; break;
 585     default: --I; --J;
 586     }
 587     for (k = 1; k <= K; ++k ) {
 588       for (j = 1; j <= J; ++j ) {
 589         for (i = 1; i <= I; ++i ) {
 590
 591           n1 = ++iNode;
 592
 593           switch (ax) {
 594           case 1: // nodes for faces normal to i direction
 595             n2 = n1 + ijLen;
 596             n3 = n2 + I;
 597             n4 = n1 + I;
 598             break;
 599           case 2: // nodes for faces normal to j direction
 600             n2 = n1 + 1;
 601             n3 = n2 + ijLen;
 602             n4 = n3 - 1;
 603             break;
 604           default: // nodes for faces normal to k direction
 605             n2 = n1 + I;
 606             n3 = n2 + 1;
 607             n4 = n1 + 1;
 608           }
 609           nodeFNumbers [ ++indexF ] = n1;
 610           nodeFNumbers [ ++indexF ] = n2;
 611           nodeFNumbers [ ++indexF ] = n3;
 612           nodeFNumbers [ ++indexF ] = n4;
 613         }
 614         if (ax != 1) ++iNode;
 615       }
 616       if (ax != 2) iNode += _iArrayLength;
 617     }
 618   }
 619   if (nbFNodes != indexF+1) {
 620     throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "Wrong nbFNodes : " \
 621                                         << nbFNodes << " indexF : " << indexF ));
 622   }
 623
 624   // fill edge node numbers
 625
 626   AX = hasEdges ? _spaceDimension : 0;
 627   for (ax = 1; ax <= AX; ++ax) {
 628
 629     iNode = 0;
 630
 631     I = _iArrayLength;
 632     J = _jArrayLength;
 633     K = _kArrayLength;
 634     if (K == 0) K = 1;
 635
 636     switch (ax) {
 637     case 1:  --I; break;
 638     case 2:  --J; break;
 639     default: --K;
 640     }
 641     for (k = 1; k <= K; ++k ) {
 642       for (j = 1; j <= J; ++j ) {
 643         for (i = 1; i <= I; ++i ) {
 644
 645           n1 = ++iNode;
 646
 647           switch (ax) {
 648           case 1: // nodes for edges going along i direction
 649             n2 = n1 + 1;
 650             break;
 651           case 2: // nodes for edges going along j direction
 652             n2 = n1 + _iArrayLength;
 653             break;
 654           default: // nodes for edges going along k direction
 655             n2 = n1 + ijLen;
 656           }
 657           nodeENumbers [ ++indexE ] = n1;
 658           nodeENumbers [ ++indexE ] = n2;
 659         }
 660         if (ax == 1) ++iNode;
 661       }
 662       if (ax == 2) iNode += _iArrayLength;
 663     }
 664   }
 665   if (nbENodes != indexE+1) {
 666     throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "Wrong nbFNodes : " \
 667                                         << nbENodes << " indexE : " << indexE ));
 668   }
 669
 670
 671   CONNECTIVITY * CellCNCT, * FaceCNCT, * EdgeCNCT;
 672
 673   // make connectivity for CELL
 674
 675   medGeometryElement aCellGeometry;
 676   if (_kArrayLength)      aCellGeometry = MED_HEXA8;
 677   else if (_jArrayLength) aCellGeometry = MED_QUAD4;
 678   else                    aCellGeometry = MED_SEG2;
 679
 680   // nodal
 681   CellCNCT = makeConnectivity (MED_CELL, aCellGeometry, nbCells, nbCNodes, nbMeshNodes, nodeCNumbers);
 682
 683   delete [] nodeCNumbers;
 684
 685   // descending
 686   {
 687     //    CellCNCT->_descending = new MEDSKYLINEARRAY() ;
 688     int * skyLineArrayIndex = new int [nbCells + 1];
 689     int nbEntitySub = CellCNCT->_type[0].getNumberOfConstituents(1);
 690     for (i=0, j=1; i <= nbCells; ++i, j += nbEntitySub)
 691       skyLineArrayIndex [ i ] = j;
 692     //    CellCNCT->_descending->setMEDSKYLINEARRAY (nbCells, nbSub,
 693     //                                               skyLineArrayIndex, subNumbers);
 694     CellCNCT->_descending = new MEDSKYLINEARRAY (nbCells, nbSub,
 695                                                  skyLineArrayIndex,
 696                                                  subNumbers);
 697     delete [] skyLineArrayIndex;
 698   }
 699   delete [] subNumbers;
 700
 701   // reverse descending
 702   {
 703     //    CellCNCT->_reverseDescendingConnectivity = new MEDSKYLINEARRAY() ;
 704     nbSub = nbRevSub/2;
 705     int * skyLineArrayIndex = new int [nbSub + 1];
 706     for (i=0, j=1; i <= nbSub; ++i, j += 2)
 707       skyLineArrayIndex [ i ] = j;
 708     //    CellCNCT->_reverseDescendingConnectivity->setMEDSKYLINEARRAY
 709     //      (nbSub, nbRevSub, skyLineArrayIndex, subRevNumbers);
 710
 711     CellCNCT->_reverseDescendingConnectivity =
 712       new MEDSKYLINEARRAY(nbSub, nbRevSub, skyLineArrayIndex, subRevNumbers);
 713
 714     delete [] skyLineArrayIndex;
 715   }
 716   delete [] subRevNumbers;
 717
 718   // make connectivity for FACE and/or EDGE
 719
 720   if (hasFaces) {
 721     FaceCNCT = makeConnectivity (MED_FACE, MED_QUAD4, nbFaces, nbFNodes, nbMeshNodes, nodeFNumbers);
 722
 723     delete [] nodeFNumbers;
 724
 725     CellCNCT->_constituent = FaceCNCT;
 726   }
 727   if (hasEdges) {
 728     EdgeCNCT = makeConnectivity (MED_EDGE, MED_SEG2, nbEdges, nbENodes, nbMeshNodes, nodeENumbers);
 729
 730     delete [] nodeENumbers;
 731
 732     if (hasFaces)
 733       FaceCNCT->_constituent = EdgeCNCT;
 734     else
 735       CellCNCT->_constituent = EdgeCNCT;
 736   }
 737
 738   MESH::_connectivity  = CellCNCT;
 739
 740   (const_cast <GRID *> (this))->_is_connectivity_filled = true;
 741
 742   END_OF(LOC);
 743 }
 744
 745 //=======================================================================
 746 //function : getArrayLength
 747 //purpose  : return array length. Axis = [1,2,3] meaning [i,j,k],
 748 //=======================================================================
 749
 750 int GRID::getArrayLength( const int Axis ) const throw (MEDEXCEPTION)
 751 {
 752   switch (Axis) {
 753   case 1: return _iArrayLength;
 754   case 2: return _jArrayLength;
 755   case 3: return _kArrayLength;
 756   default:
 757     throw MED_EXCEPTION ( LOCALIZED( STRING("GRID::getArrayLength ( ") << Axis << ")"));
 758   }
 759   return 0;
 760 }
 761
 762 //=======================================================================
 763 //function : getArrayValue
 764 //purpose  : return i-th array component. Axis = [1,2,3] meaning [i,j,k],
 765 //           exception if Axis out of [1-3] range
 766 //           exception if i is out of range 0 <= i < getArrayLength(Axis);
 767 //=======================================================================
 768
 769 const double GRID::getArrayValue (const int Axis, const int i) const throw (MEDEXCEPTION)
 770 {
 771   if (i < 0 || i >= getArrayLength(Axis))
 772     throw MED_EXCEPTION
 773       ( LOCALIZED(STRING("GRID::getArrayValue ( ") << Axis << ", " << i << ")"));
 774
 775   switch (Axis) {
 776   case 1:  return _iArray[ i ];
 777   case 2:  return _jArray[ i ];
 778   default: return _kArray[ i ];
 779   }
 780   return 0.0;
 781 }
 782
 783 //=======================================================================
 784 //function : getEdgeNumber
 785 //purpose  :
 786 //=======================================================================
 787
 788 int GRID::getEdgeNumber(const int Axis, const int i, const int j, const int k)
 789   const throw (MEDEXCEPTION)
 790 {
 791   const char * LOC = "GRID::getEdgeNumber(Axis, i,j,k) :";
 792
 793   BEGIN_OF(LOC);
 794
 795   int Len[4] = {0,_iArrayLength, _jArrayLength, _kArrayLength }, I=1, J=2, K=3;
 796   int maxAxis = Len[ K ] ? 3 : 2;
 797
 798   if (Axis <= 0 || Axis > maxAxis)
 799     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Axis out of range: " << Axis));
 800
 801   Len[Axis]--;
 802   int Nb = 1 + i + j*Len[ I ] + k*Len[ J ]*Len[ K ];
 803   Len[Axis]++ ;
 804
 805   if (!Len[ K ])
 806     Len[ K ] = 1;
 807
 808   if (Axis > 1) { // add all edges in i direction
 809     Len[I]-- ;
 810     Nb += Len[ I ]*Len[ J ]*Len[ K ];
 811     Len[I]++ ;
 812   }
 813   if (Axis > 2) { // add all edges in j direction
 814     Len[J]-- ;
 815     Nb += Len[ I ]*Len[ J ]*Len[ K ];
 816   }
 817
 818   END_OF(LOC);
 819
 820   return Nb;
 821 }
 822
 823 //=======================================================================
 824 //function : getFaceNumber
 825 //purpose  : return a NODE, EDGE, FACE, CELL number by its position in the grid.
 826 //           Axis [1,2,3] means one of directions: along i, j or k
 827 //           For Cell contituents (FACE or EDGE), Axis selects one of those having same (i,j,k):
 828 //           * a FACE which is normal to direction along given Axis;
 829 //           * an EDGE going along given Axis.
 830 //           Exception for Axis out of range
 831 //=======================================================================
 832
 833 int GRID::getFaceNumber(const int Axis, const int i, const int j, const int k)
 834   const throw (MEDEXCEPTION)
 835 {
 836   const char * LOC = "GRID::getFaceNumber(Axis, i,j,k) :";
 837
 838   BEGIN_OF(LOC);
 839
 840 //  if (Axis <= 0 || Axis > 3)
 841   if (Axis < 0 || Axis > 3)
 842     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Axis = " << Axis));
 843
 844   int Len[4] = {0,_iArrayLength-1, _jArrayLength-1, _kArrayLength-1 }, I=1, J=2, K=3;
 845
 846   Len[Axis]++;
 847   int Nb = 1 + i + j*Len[ I ] + k*Len[ J ]*Len[ K ];
 848   Len[Axis]--;
 849
 850   if (Axis > 1) { // add all faces in i direction
 851     Len[I]++ ;
 852     Nb += Len[ I ]*Len[ J ]*Len[ K ];
 853     Len[I]-- ;
 854   }
 855   if (Axis > 2) { // add all faces in j direction
 856     Len[J]++ ;
 857     Nb += Len[ I ]*Len[ J ]*Len[ K ];
 858   }
 859
 860   END_OF(LOC);
 861
 862   return Nb;
 863 }
 864
 865 //=======================================================================
 866 //function : getNodePosition
 867 //purpose  :
 868 //=======================================================================
 869
 870 void GRID::getNodePosition(const int Number, int& i, int& j, int& k) const
 871   throw (MEDEXCEPTION)
 872 {
 873   const char * LOC = "GRID::getNodePosition(Number, i,j,k) :";
 874
 875   BEGIN_OF(LOC);
 876
 877   if (Number <= 0 || Number > _numberOfNodes)
 878     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
 879
 880   int Len[] = {_iArrayLength, _jArrayLength, _kArrayLength }, I=0, J=1;
 881   // , K=2; !! UNUSED VARIABLE !!
 882
 883   int ijLen = Len[I] * Len[J]; // nb in a full k layer
 884   int kLen = (Number - 1) % ijLen;    // nb in the non full k layer
 885
 886   i = kLen % Len[J];
 887   j = kLen / Len[J];
 888   k = (Number - 1) / ijLen;
 889
 890   ////cout <<" NODE POS: " << Number << " - " << i << ", " << j << ", " << k << endl;
 891
 892   END_OF(LOC);
 893
 894 }
 895
 896 //=======================================================================
 897 //function : getCellPosition
 898 //purpose  :
 899 //=======================================================================
 900
 901 void GRID::getCellPosition(const int Number, int& i, int& j, int& k) const
 902   throw (MEDEXCEPTION)
 903 {
 904   const char * LOC = "GRID::getCellPosition(Number, i,j,k) :";
 905
 906   BEGIN_OF(LOC);
 907
 908   int Len[4] = {0,_iArrayLength-1, _jArrayLength-1, _kArrayLength-1 }, I=1, J=2;
 909   // , K=3; !! UNUSED VARIABLE !!
 910
 911 //  if (Number <= 0 || Number > getCellNumber(Len[I]-1, Len[J]-1, Len[K]-1))
 912 //    throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
 913
 914   int ijLen = Len[I] * Len[J]; // nb in a full k layer
 915   int kLen = (Number - 1) % ijLen;    // nb in the non full k layer
 916
 917   i = kLen % Len[J];
 918   j = kLen / Len[J];
 919   k = (Number - 1) / ijLen;
 920
 921   END_OF(LOC);
 922 }
 923
 924 //=======================================================================
 925 //function : getEdgePosition
 926 //purpose  :
 927 //=======================================================================
 928
 929 void GRID::getEdgePosition(const int Number, int& Axis, int& i, int& j, int& k)
 930   const throw (MEDEXCEPTION)
 931 {
 932   const char * LOC = "GRID::getEdgePosition(Number, i,j,k) :";
 933
 934   BEGIN_OF(LOC);
 935
 936   if (!_jArrayLength)
 937     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "no edges in the grid: "));
 938
 939   if (Number <= 0)
 940     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
 941
 942   int Len[4] = {0,_iArrayLength, _jArrayLength, _kArrayLength }, I=1, J=2, K=3;
 943
 944   int theNb = Number;
 945   int maxAxis = _kArrayLength ? 3 : 2;
 946
 947   for (Axis = 1; Axis <= maxAxis; ++Axis)
 948   {
 949     Len[Axis]--;  // one less edge in Axis direction
 950
 951     // max edge number in Axis direction
 952     int maxNb = getEdgeNumber (Axis, Len[I]-1, Len[J]-1, Len[K]-1);
 953
 954     if (theNb > maxNb)
 955     {
 956       Len[Axis]++;
 957       theNb -= maxNb;
 958       continue;
 959     }
 960
 961     if (theNb == maxNb)
 962     {
 963       i = Len[I]-1;
 964       j = Len[J]-1;
 965       k = Len[K]-1;
 966     }
 967     else
 968     {
 969       int ijLen = Len[I] * Len[J]; // nb in a full k layer
 970       int kLen = (theNb - 1) % ijLen;    // nb in the non full k layer
 971       i = kLen % Len[J];
 972       j = kLen / Len[J];
 973       k = (theNb - 1) / ijLen;
 974     }
 975
 976     END_OF(LOC);
 977
 978     return;
 979   }
 980
 981   throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
 982 }
 983
 984 //=======================================================================
 985 //function : getFacePosition
 986 //purpose  : return position (i,j,k) of an entity #Number
 987 //           Axis [1,2,3] means one of directions: along i, j or k
 988 //           For Cell contituents (FACE or EDGE), Axis selects one of those having same (i,j,k):
 989 //           * a FACE which is normal to direction along given Axis;
 990 //           * an EDGE going along given Axis.
 991 //           Exception for Number out of range
 992 //=======================================================================
 993
 994 void GRID::getFacePosition(const int Number, int& Axis, int& i, int& j, int& k)
 995   const throw (MEDEXCEPTION)
 996 {
 997   const char * LOC = "GRID::getFacePosition(Number, i,j,k) :";
 998
 999   BEGIN_OF(LOC);
1000
1001   if (_kArrayLength == 0) {
1002     getCellPosition(Number, i, j, k);
1003     Axis = 1;
1004     return;
1005   };
1006
1007   if (!_kArrayLength)
1008     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "no faces in the grid: "));
1009
1010   if ( Number <= 0 )
1011     throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
1012
1013   int Len[4] = {0,_iArrayLength-1, _jArrayLength-1, _kArrayLength-1 }, I=1, J=2, K=3;
1014   int theNb = Number;
1015
1016   for (Axis = 1; Axis <= 3; ++Axis)
1017   {
1018     Len[Axis]++;
1019
1020     // max face number in Axis direction
1021     int maxNb = getFaceNumber (Axis, Len[I]-1, Len[J]-1, Len[K]-1);
1022
1023     if (theNb > maxNb)
1024     {
1025       Len[Axis]--;
1026       theNb -= maxNb;
1027       continue;
1028     }
1029
1030     if (theNb == maxNb)
1031     {
1032       i = Len[I]-1;
1033       j = Len[J]-1;
1034       k = Len[K]-1;
1035     }
1036     else
1037     {
1038       int ijLen = Len[I] * Len[J]; // nb in a full k layer
1039       int kLen = (theNb - 1)  % ijLen;    // nb in the non full k layer
1040       i = kLen % Len[J];
1041       j = kLen / Len[J];
1042       k = (theNb - 1) / ijLen;
1043     }
1044
1045     END_OF(LOC);
1046
1047     return;
1048   }
1049
1050   throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
1051 }
1052
1053 //=======================================================================
1054 //function : writeUnstructured
1055 //purpose  : write a Grid as an Unstructured mesh
1056 //=======================================================================
1057
1058 void GRID::writeUnstructured(int index, const string & driverName)
1059 {
1060   const char * LOC = "GRID::writeUnstructured(int index=0, const string & driverName) : ";
1061   BEGIN_OF(LOC);
1062
1063   if ( _drivers[index] ) {
1064
1065     makeUnstructured();
1066     _isAGrid = false;
1067
1068     _drivers[index]->open();
1069     if (driverName != "") _drivers[index]->setMeshName(driverName);
1070     _drivers[index]->write();
1071     _drivers[index]->close();
1072
1073     _isAGrid = true;
1074   }
1075   else
1076     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
1077                                     << "The index given is invalid, index must be between  0 and |"
1078                                     << _drivers.size()
1079                                     )
1080                          );
1081   END_OF(LOC);
1082 }
1083
1084 void GRID::read(int index)
1085 {
1086   const char * LOC = "GRID::read(int index=0) : ";
1087   BEGIN_OF(LOC);
1088
1089   if (_drivers[index]) {
1090     _drivers[index]->open();
1091     _drivers[index]->read();
1092     _drivers[index]->close();
1093   }
1094   else
1095     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
1096                                      << "The index given is invalid, index must be between  0 and |"
1097                                      << _drivers.size()
1098                                      )
1099                           );
1100   if (_isAGrid)
1101     fillMeshAfterRead();
1102
1103   END_OF(LOC);
1104 }