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