src/MEDMEM/MEDMEM_Grid.cxx

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