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