src/MEDMEM/MEDMEM_Grid.cxx

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