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