src/MEDMEM/MEDMEM_Grid.cxx

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