src/MEDMEM/MEDMEM_Grid.cxx

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