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

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