applications/utilities/mesh/conversion/plot3dToFoam/hexBlock.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "hexBlock.H"
  28
  29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33
  34 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  35
  36 label hexBlock::vtxLabel(label a, label b, label c) const
  37 {
  38     return (a + b*(xDim_ + 1) + c*(xDim_ + 1)*(yDim_ + 1));
  39 }
  40
  41
  42 // Calculate the handedness of the block by looking at the orientation
  43 // of the spanning edges of a hex. Loops until valid cell found (since might
  44 // be prism)
  45 void hexBlock::setHandedness()
  46 {
  47     const pointField& p = points_;
  48
  49     for (label k = 0; k <= zDim_ - 1; k++)
  50     {
  51         for (label j = 0; j <= yDim_ - 1; j++)
  52         {
  53             for (label i = 0; i <= xDim_ - 1; i++)
  54             {
  55                 vector x = p[vtxLabel(i+1, j, k)] - p[vtxLabel(i, j, k)];
  56                 vector y = p[vtxLabel(i, j+1, k)] - p[vtxLabel(i, j, k)];
  57                 vector z = p[vtxLabel(i, j, k+1)] - p[vtxLabel(i, j, k)];
  58
  59                 if (mag(x) > SMALL && mag(y) > SMALL && mag(z) > SMALL)
  60                 {
  61                     Info<< "Looking at cell "
  62                         << i << ' ' << j << ' ' << k
  63                         << " to determine orientation."
  64                         << endl;
  65
  66                     if (((x ^ y) & z) > 0)
  67                     {
  68                         Info<< "Right-handed block." << endl;
  69                         blockHandedness_ = right;
  70                     }
  71                     else
  72                     {
  73                         Info << "Left-handed block." << endl;
  74                         blockHandedness_ = left;
  75                     }
  76                     return;
  77                 }
  78                 else
  79                 {
  80                     Info<< "Cannot determine orientation of cell "
  81                         << i << ' ' << j << ' ' << k
  82                         << " since has base vectors " << x << y << z << endl;
  83                 }
  84             }
  85         }
  86     }
  87
  88     if (blockHandedness_ == noPoints)
  89     {
  90         WarningIn("hexBlock::hexBlock::setHandedness()")
  91             << "Cannot determine orientation of block."
  92             << " Continuing as if right handed." << endl;
  93         blockHandedness_ = right;
  94     }
  95 }
  96
  97
  98 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  99
 100 // Construct from components
 101 hexBlock::hexBlock(const label nx, const label ny, const label nz)
 102 :
 103     xDim_(nx - 1),
 104     yDim_(ny - 1),
 105     zDim_(nz - 1),
 106     blockHandedness_(noPoints),
 107     points_((xDim_ + 1)*(yDim_ + 1)*(zDim_ + 1))
 108 {}
 109
 110
 111 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 112
 113 void hexBlock::readPoints
 114 (
 115     const bool readBlank,
 116     const scalar twoDThicknes,
 117     Istream& is
 118 )
 119 {
 120     scalar iBlank;
 121
 122     label nPoints = points_.size();
 123
 124     if (twoDThicknes > 0)
 125     {
 126         nPoints /= 2;
 127     }
 128
 129     Info<< "Reading " << nPoints << " x coordinates..." << endl;
 130     for (label i=0; i < nPoints; i++)
 131     {
 132         is >> points_[i].x();
 133     }
 134
 135     Info<< "Reading " << nPoints << " y coordinates..." << endl;
 136     for (label i=0; i < nPoints; i++)
 137     {
 138         is >> points_[i].y();
 139     }
 140
 141     if (twoDThicknes > 0)
 142     {
 143         Info<< "Extruding " << nPoints << " points in z direction..." << endl;
 144         // Duplicate points
 145         for (label i=0; i < nPoints; i++)
 146         {
 147             points_[i+nPoints] = points_[i];
 148         }
 149         for (label i=0; i < nPoints; i++)
 150         {
 151             points_[i].z() = 0;
 152             points_[i+nPoints].z() = twoDThicknes;
 153         }
 154     }
 155     else
 156     {
 157         Info<< "Reading " << nPoints << " z coordinates..." << endl;
 158         for (label i=0; i < nPoints; i++)
 159         {
 160             is >> points_[i].z();
 161         }
 162     }
 163
 164
 165     if (readBlank)
 166     {
 167         Info<< "Reading " << nPoints << " blanks..." << endl;
 168         for (label i=0; i < nPoints; i++)
 169         {
 170             is >> iBlank;
 171         }
 172     }
 173
 174     // Set left- or righthandedness of block by looking at a cell.
 175     setHandedness();
 176 }
 177
 178
 179 labelListList hexBlock::blockCells() const
 180 {
 181     labelListList result(xDim_*yDim_*zDim_);
 182
 183     label cellNo = 0;
 184
 185     if (blockHandedness_ == right)
 186     {
 187         for (label k = 0; k <= zDim_ - 1; k++)
 188         {
 189             for (label j = 0; j <= yDim_ - 1; j++)
 190             {
 191                 for (label i = 0; i <= xDim_ - 1; i++)
 192                 {
 193                     labelList& hexLabels = result[cellNo];
 194                     hexLabels.setSize(8);
 195
 196                     hexLabels[0] = vtxLabel(i, j, k);
 197                     hexLabels[1] = vtxLabel(i+1, j, k);
 198                     hexLabels[2] = vtxLabel(i+1, j+1, k);
 199                     hexLabels[3] = vtxLabel(i, j+1, k);
 200                     hexLabels[4] = vtxLabel(i, j, k+1);
 201                     hexLabels[5] = vtxLabel(i+1, j, k+1);
 202                     hexLabels[6] = vtxLabel(i+1, j+1, k+1);
 203                     hexLabels[7] = vtxLabel(i, j+1, k+1);
 204
 205                     cellNo++;
 206                 }
 207             }
 208         }
 209     }
 210     else if (blockHandedness_ == left)
 211     {
 212         for (label k = 0; k <= zDim_ - 1; k++)
 213         {
 214             for (label j = 0; j <= yDim_ - 1; j++)
 215             {
 216                 for (label i = 0; i <= xDim_ - 1; i++)
 217                 {
 218                     labelList& hexLabels = result[cellNo];
 219                     hexLabels.setSize(8);
 220
 221                     hexLabels[0] = vtxLabel(i, j, k+1);
 222                     hexLabels[1] = vtxLabel(i+1, j, k+1);
 223                     hexLabels[2] = vtxLabel(i+1, j+1, k+1);
 224                     hexLabels[3] = vtxLabel(i, j+1, k+1);
 225                     hexLabels[4] = vtxLabel(i, j, k);
 226                     hexLabels[5] = vtxLabel(i+1, j, k);
 227                     hexLabels[6] = vtxLabel(i+1, j+1, k);
 228                     hexLabels[7] = vtxLabel(i, j+1, k);
 229
 230                     cellNo++;
 231                 }
 232             }
 233         }
 234     }
 235     else
 236     {
 237         FatalErrorIn("hexBlock::cellShapes()")
 238             << "Unable to determine block handedness as points "
 239             << "have not been read in yet"
 240             << abort(FatalError);
 241     }
 242
 243     return result;
 244 }
 245
 246
 247 // Return block patch faces given direction and range limits
 248 // From the cfx manual: direction
 249 // 0 = solid (3-D patch),
 250 // 1 = high i, 2 = high j, 3 = high k
 251 // 4 = low i, 5 = low j, 6 = low k
 252 faceList hexBlock::patchFaces(const label direc, const labelList& range) const
 253 {
 254     if (range.size() != 6)
 255     {
 256         FatalErrorIn
 257         (
 258             "patchFaces(const label direc, const labelList& range) const"
 259         )   << "Invalid size of the range array: " << range.size()
 260             << ". Should be 6 (xMin, xMax, yMin, yMax, zMin, zMax"
 261             << abort(FatalError);
 262     }
 263
 264     label xMinRange = range[0];
 265     label xMaxRange = range[1];
 266     label yMinRange = range[2];
 267     label yMaxRange = range[3];
 268     label zMinRange = range[4];
 269     label zMaxRange = range[5];
 270
 271     faceList result(0);
 272
 273     switch (direc)
 274     {
 275         case 1:
 276         {
 277             // high i = xmax
 278
 279             result.setSize
 280             (
 281                 (yMaxRange - yMinRange + 1)*(zMaxRange - zMinRange + 1)
 282             );
 283
 284             label p = 0;
 285             for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
 286             {
 287                 for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
 288                 {
 289                     result[p].setSize(4);
 290
 291                     // set the points
 292                     result[p][0] = vtxLabel(xDim_, j, k);
 293                     result[p][1] = vtxLabel(xDim_, j+1, k);
 294                     result[p][2] = vtxLabel(xDim_, j+1, k+1);
 295                     result[p][3] = vtxLabel(xDim_, j, k+1);
 296
 297                     p++;
 298                 }
 299             }
 300
 301             result.setSize(p);
 302             break;
 303         }
 304
 305         case 2:
 306         {
 307             // high j = ymax
 308             result.setSize
 309             (
 310                 (xMaxRange - xMinRange + 1)*(zMaxRange - zMinRange + 1)
 311             );
 312
 313             label p = 0;
 314             for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
 315             {
 316                 for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
 317                 {
 318                     result[p].setSize(4);
 319
 320                     // set the points
 321                     result[p][0] = vtxLabel(i, yDim_, k);
 322                     result[p][1] = vtxLabel(i, yDim_, k + 1);
 323                     result[p][2] = vtxLabel(i + 1, yDim_, k + 1);
 324                     result[p][3] = vtxLabel(i + 1, yDim_, k);
 325
 326                     p++;
 327                 }
 328             }
 329
 330             result.setSize(p);
 331             break;
 332         }
 333
 334         case 3:
 335         {
 336             // high k = zmax
 337             result.setSize
 338             (
 339                 (xMaxRange - xMinRange + 1)*(yMaxRange - yMinRange + 1)
 340             );
 341
 342             label p = 0;
 343             for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
 344             {
 345                 for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
 346                 {
 347                     result[p].setSize(4);
 348
 349                     // set the points
 350                     result[p][0] = vtxLabel(i, j, zDim_);
 351                     result[p][1] = vtxLabel(i + 1, j, zDim_);
 352                     result[p][2] = vtxLabel(i + 1, j + 1, zDim_);
 353                     result[p][3] = vtxLabel(i, j + 1, zDim_);
 354
 355                     p++;
 356                 }
 357             }
 358
 359             result.setSize(p);
 360             break;
 361         }
 362
 363         case 4:
 364         {
 365             // low i = xmin
 366             result.setSize
 367             (
 368                 (yMaxRange - yMinRange + 1)*(zMaxRange - zMinRange + 1)
 369             );
 370
 371             label p = 0;
 372             for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
 373             {
 374                 for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
 375                 {
 376                     result[p].setSize(4);
 377
 378                     // set the points
 379                     result[p][0] = vtxLabel(0, j, k);
 380                     result[p][1] = vtxLabel(0, j, k + 1);
 381                     result[p][2] = vtxLabel(0, j + 1, k + 1);
 382                     result[p][3] = vtxLabel(0, j + 1, k);
 383
 384                     p++;
 385                 }
 386             }
 387
 388             result.setSize(p);
 389             break;
 390         }
 391
 392         case 5:
 393         {
 394             // low j = ymin
 395             result.setSize
 396             (
 397                 (xMaxRange - xMinRange + 1)*(zMaxRange - zMinRange + 1)
 398             );
 399
 400             label p = 0;
 401             for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
 402             {
 403                 for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
 404                 {
 405                     result[p].setSize(4);
 406
 407                     // set the points
 408                     result[p][0] = vtxLabel(i, 0, k);
 409                     result[p][1] = vtxLabel(i + 1, 0, k);
 410                     result[p][2] = vtxLabel(i + 1, 0, k + 1);
 411                     result[p][3] = vtxLabel(i, 0, k + 1);
 412
 413                     p++;
 414                 }
 415             }
 416
 417             result.setSize(p);
 418             break;
 419         }
 420
 421         case 6:
 422         {
 423             // low k = zmin
 424             result.setSize
 425             (
 426                 (xMaxRange - xMinRange + 1)*(yMaxRange - yMinRange + 1)
 427             );
 428
 429             label p = 0;
 430             for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
 431             {
 432                 for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
 433                 {
 434                     result[p].setSize(4);
 435
 436                     // set the points
 437                     result[p][0] = vtxLabel(i, j, 0);
 438                     result[p][1] = vtxLabel(i, j + 1, 0);
 439                     result[p][2] = vtxLabel(i + 1, j + 1, 0);
 440                     result[p][3] = vtxLabel(i + 1, j, 0);
 441
 442                     p++;
 443                 }
 444             }
 445
 446             result.setSize(p);
 447             break;
 448         }
 449
 450         default:
 451         {
 452             FatalErrorIn
 453             (
 454                 "patchFaces(const label direc, const labelList& range) const"
 455             )   << "direction out of range (1 to 6): " << direc
 456                 << abort(FatalError);
 457         }
 458     }
 459
 460     // Correct the face orientation based on the handedness of the block.
 461     // Do nothing for the right-handed block
 462     if (blockHandedness_ == noPoints)
 463     {
 464         FatalErrorIn
 465         (
 466             "patchFaces(const label direc, const labelList& range) const"
 467         )   << "Unable to determine block handedness as points "
 468             << "have not been read in yet"
 469             << abort(FatalError);
 470     }
 471     else if (blockHandedness_ == left)
 472     {
 473         // turn all faces inside out
 474         forAll (result, faceI)
 475         {
 476             result[faceI] = result[faceI].reverseFace();
 477         }
 478     }
 479
 480     return result;
 481 }
 482
 483
 484 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 485
 486 } // End namespace Foam
 487
 488 // ************************************************************************* //