FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

hexBlock.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "hexBlock.H"
00027 
00028 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00029 
00030 namespace Foam
00031 {
00032 
00033 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00034 
00035 label hexBlock::vtxLabel(label a, label b, label c) const
00036 {
00037     return (a + b*(xDim_ + 1) + c*(xDim_ + 1)*(yDim_ + 1));
00038 }
00039 
00040 
00041 // Calculate the handedness of the block by looking at the orientation
00042 // of the spanning edges of a hex. Loops until valid cell found (since might
00043 // be prism)
00044 void hexBlock::setHandedness()
00045 {
00046     const pointField& p = points_;
00047 
00048     for (label k = 0; k <= zDim_ - 1; k++)
00049     {
00050         for (label j = 0; j <= yDim_ - 1; j++)
00051         {
00052             for (label i = 0; i <= xDim_ - 1; i++)
00053             {
00054                 vector x = p[vtxLabel(i+1, j, k)] - p[vtxLabel(i, j, k)];
00055                 vector y = p[vtxLabel(i, j+1, k)] - p[vtxLabel(i, j, k)];
00056                 vector z = p[vtxLabel(i, j, k+1)] - p[vtxLabel(i, j, k)];
00057 
00058                 if (mag(x) > SMALL && mag(y) > SMALL && mag(z) > SMALL)
00059                 {
00060                     Info<< "Looking at cell "
00061                         << i << ' ' << j << ' ' << k
00062                         << " to determine orientation."
00063                         << endl;
00064 
00065                     if (((x ^ y) & z) > 0)
00066                     {
00067                         Info<< "Right-handed block." << endl;
00068                         blockHandedness_ = right;
00069                     }
00070                     else
00071                     {
00072                         Info << "Left-handed block." << endl;
00073                         blockHandedness_ = left;
00074                     }
00075                     return;
00076                 }
00077                 else
00078                 {
00079                     Info<< "Cannot determine orientation of cell "
00080                         << i << ' ' << j << ' ' << k
00081                         << " since has base vectors " << x << y << z << endl;
00082                 }
00083             }
00084         }
00085     }
00086 
00087     if (blockHandedness_ == noPoints)
00088     {
00089         WarningIn("hexBlock::hexBlock::setHandedness()")
00090             << "Cannot determine orientation of block."
00091             << " Continuing as if right handed." << endl;
00092         blockHandedness_ = right;
00093     }
00094 }
00095 
00096 
00097 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00098 
00099 // Construct from components
00100 hexBlock::hexBlock(const label nx, const label ny, const label nz)
00101 :
00102     xDim_(nx - 1),
00103     yDim_(ny - 1),
00104     zDim_(nz - 1),
00105     blockHandedness_(noPoints),
00106     points_((xDim_ + 1)*(yDim_ + 1)*(zDim_ + 1))
00107 {}
00108 
00109 
00110 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00111 
00112 void hexBlock::readPoints
00113 (
00114     const bool readBlank,
00115     const scalar twoDThicknes,
00116     Istream& is
00117 )
00118 {
00119     scalar iBlank;
00120 
00121     label nPoints = points_.size();
00122 
00123     if (twoDThicknes > 0)
00124     {
00125         nPoints /= 2;
00126     }
00127 
00128     Info<< "Reading " << nPoints << " x coordinates..." << endl;
00129     for (label i=0; i < nPoints; i++)
00130     {
00131         is >> points_[i].x();
00132     }
00133 
00134     Info<< "Reading " << nPoints << " y coordinates..." << endl;
00135     for (label i=0; i < nPoints; i++)
00136     {
00137         is >> points_[i].y();
00138     }
00139 
00140     if (twoDThicknes > 0)
00141     {
00142         Info<< "Extruding " << nPoints << " points in z direction..." << endl;
00143         // Duplicate points
00144         for (label i=0; i < nPoints; i++)
00145         {
00146             points_[i+nPoints] = points_[i];
00147         }
00148         for (label i=0; i < nPoints; i++)
00149         {
00150             points_[i].z() = 0;
00151             points_[i+nPoints].z() = twoDThicknes;
00152         }
00153     }
00154     else
00155     {
00156         Info<< "Reading " << nPoints << " z coordinates..." << endl;
00157         for (label i=0; i < nPoints; i++)
00158         {
00159             is >> points_[i].z();
00160         }
00161     }
00162 
00163 
00164     if (readBlank)
00165     {
00166         Info<< "Reading " << nPoints << " blanks..." << endl;
00167         for (label i=0; i < nPoints; i++)
00168         {
00169             is >> iBlank;
00170         }
00171     }
00172 
00173     // Set left- or righthandedness of block by looking at a cell.
00174     setHandedness();
00175 }
00176 
00177 
00178 labelListList hexBlock::blockCells() const
00179 {
00180     labelListList result(xDim_*yDim_*zDim_);
00181 
00182     label cellNo = 0;
00183 
00184     if (blockHandedness_ == right)
00185     {
00186         for (label k = 0; k <= zDim_ - 1; k++)
00187         {
00188             for (label j = 0; j <= yDim_ - 1; j++)
00189             {
00190                 for (label i = 0; i <= xDim_ - 1; i++)
00191                 {
00192                     labelList& hexLabels = result[cellNo];
00193                     hexLabels.setSize(8);
00194 
00195                     hexLabels[0] = vtxLabel(i, j, k);
00196                     hexLabels[1] = vtxLabel(i+1, j, k);
00197                     hexLabels[2] = vtxLabel(i+1, j+1, k);
00198                     hexLabels[3] = vtxLabel(i, j+1, k);
00199                     hexLabels[4] = vtxLabel(i, j, k+1);
00200                     hexLabels[5] = vtxLabel(i+1, j, k+1);
00201                     hexLabels[6] = vtxLabel(i+1, j+1, k+1);
00202                     hexLabels[7] = vtxLabel(i, j+1, k+1);
00203 
00204                     cellNo++;
00205                 }
00206             }
00207         }
00208     }
00209     else if (blockHandedness_ == left)
00210     {
00211         for (label k = 0; k <= zDim_ - 1; k++)
00212         {
00213             for (label j = 0; j <= yDim_ - 1; j++)
00214             {
00215                 for (label i = 0; i <= xDim_ - 1; i++)
00216                 {
00217                     labelList& hexLabels = result[cellNo];
00218                     hexLabels.setSize(8);
00219 
00220                     hexLabels[0] = vtxLabel(i, j, k+1);
00221                     hexLabels[1] = vtxLabel(i+1, j, k+1);
00222                     hexLabels[2] = vtxLabel(i+1, j+1, k+1);
00223                     hexLabels[3] = vtxLabel(i, j+1, k+1);
00224                     hexLabels[4] = vtxLabel(i, j, k);
00225                     hexLabels[5] = vtxLabel(i+1, j, k);
00226                     hexLabels[6] = vtxLabel(i+1, j+1, k);
00227                     hexLabels[7] = vtxLabel(i, j+1, k);
00228 
00229                     cellNo++;
00230                 }
00231             }
00232         }
00233     }
00234     else
00235     {
00236         FatalErrorIn("hexBlock::cellShapes()")
00237             << "Unable to determine block handedness as points "
00238             << "have not been read in yet"
00239             << abort(FatalError);
00240     }
00241 
00242     return result;
00243 }
00244 
00245 
00246 // Return block patch faces given direction and range limits
00247 // From the cfx manual: direction
00248 // 0 = solid (3-D patch),
00249 // 1 = high i, 2 = high j, 3 = high k
00250 // 4 = low i, 5 = low j, 6 = low k
00251 faceList hexBlock::patchFaces(const label direc, const labelList& range) const
00252 {
00253     if (range.size() != 6)
00254     {
00255         FatalErrorIn
00256         (
00257             "patchFaces(const label direc, const labelList& range) const"
00258         )   << "Invalid size of the range array: " << range.size()
00259             << ". Should be 6 (xMin, xMax, yMin, yMax, zMin, zMax"
00260             << abort(FatalError);
00261     }
00262 
00263     label xMinRange = range[0];
00264     label xMaxRange = range[1];
00265     label yMinRange = range[2];
00266     label yMaxRange = range[3];
00267     label zMinRange = range[4];
00268     label zMaxRange = range[5];
00269 
00270     faceList result(0);
00271 
00272     switch (direc)
00273     {
00274         case 1:
00275         {
00276             // high i = xmax
00277 
00278             result.setSize
00279             (
00280                 (yMaxRange - yMinRange + 1)*(zMaxRange - zMinRange + 1)
00281             );
00282 
00283             label p = 0;
00284             for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
00285             {
00286                 for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
00287                 {
00288                     result[p].setSize(4);
00289 
00290                     // set the points
00291                     result[p][0] = vtxLabel(xDim_, j, k);
00292                     result[p][1] = vtxLabel(xDim_, j+1, k);
00293                     result[p][2] = vtxLabel(xDim_, j+1, k+1);
00294                     result[p][3] = vtxLabel(xDim_, j, k+1);
00295 
00296                     p++;
00297                 }
00298             }
00299 
00300             result.setSize(p);
00301             break;
00302         }
00303 
00304         case 2:
00305         {
00306             // high j = ymax
00307             result.setSize
00308             (
00309                 (xMaxRange - xMinRange + 1)*(zMaxRange - zMinRange + 1)
00310             );
00311 
00312             label p = 0;
00313             for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
00314             {
00315                 for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
00316                 {
00317                     result[p].setSize(4);
00318 
00319                     // set the points
00320                     result[p][0] = vtxLabel(i, yDim_, k);
00321                     result[p][1] = vtxLabel(i, yDim_, k + 1);
00322                     result[p][2] = vtxLabel(i + 1, yDim_, k + 1);
00323                     result[p][3] = vtxLabel(i + 1, yDim_, k);
00324 
00325                     p++;
00326                 }
00327             }
00328 
00329             result.setSize(p);
00330             break;
00331         }
00332 
00333         case 3:
00334         {
00335             // high k = zmax
00336             result.setSize
00337             (
00338                 (xMaxRange - xMinRange + 1)*(yMaxRange - yMinRange + 1)
00339             );
00340 
00341             label p = 0;
00342             for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
00343             {
00344                 for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
00345                 {
00346                     result[p].setSize(4);
00347 
00348                     // set the points
00349                     result[p][0] = vtxLabel(i, j, zDim_);
00350                     result[p][1] = vtxLabel(i + 1, j, zDim_);
00351                     result[p][2] = vtxLabel(i + 1, j + 1, zDim_);
00352                     result[p][3] = vtxLabel(i, j + 1, zDim_);
00353 
00354                     p++;
00355                 }
00356             }
00357 
00358             result.setSize(p);
00359             break;
00360         }
00361 
00362         case 4:
00363         {
00364             // low i = xmin
00365             result.setSize
00366             (
00367                 (yMaxRange - yMinRange + 1)*(zMaxRange - zMinRange + 1)
00368             );
00369 
00370             label p = 0;
00371             for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
00372             {
00373                 for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
00374                 {
00375                     result[p].setSize(4);
00376 
00377                     // set the points
00378                     result[p][0] = vtxLabel(0, j, k);
00379                     result[p][1] = vtxLabel(0, j, k + 1);
00380                     result[p][2] = vtxLabel(0, j + 1, k + 1);
00381                     result[p][3] = vtxLabel(0, j + 1, k);
00382 
00383                     p++;
00384                 }
00385             }
00386 
00387             result.setSize(p);
00388             break;
00389         }
00390 
00391         case 5:
00392         {
00393             // low j = ymin
00394             result.setSize
00395             (
00396                 (xMaxRange - xMinRange + 1)*(zMaxRange - zMinRange + 1)
00397             );
00398 
00399             label p = 0;
00400             for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
00401             {
00402                 for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
00403                 {
00404                     result[p].setSize(4);
00405 
00406                     // set the points
00407                     result[p][0] = vtxLabel(i, 0, k);
00408                     result[p][1] = vtxLabel(i + 1, 0, k);
00409                     result[p][2] = vtxLabel(i + 1, 0, k + 1);
00410                     result[p][3] = vtxLabel(i, 0, k + 1);
00411 
00412                     p++;
00413                 }
00414             }
00415 
00416             result.setSize(p);
00417             break;
00418         }
00419 
00420         case 6:
00421         {
00422             // low k = zmin
00423             result.setSize
00424             (
00425                 (xMaxRange - xMinRange + 1)*(yMaxRange - yMinRange + 1)
00426             );
00427 
00428             label p = 0;
00429             for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
00430             {
00431                 for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
00432                 {
00433                     result[p].setSize(4);
00434 
00435                     // set the points
00436                     result[p][0] = vtxLabel(i, j, 0);
00437                     result[p][1] = vtxLabel(i, j + 1, 0);
00438                     result[p][2] = vtxLabel(i + 1, j + 1, 0);
00439                     result[p][3] = vtxLabel(i + 1, j, 0);
00440 
00441                     p++;
00442                 }
00443             }
00444 
00445             result.setSize(p);
00446             break;
00447         }
00448 
00449         default:
00450         {
00451             FatalErrorIn
00452             (
00453                 "patchFaces(const label direc, const labelList& range) const"
00454             )   << "direction out of range (1 to 6): " << direc
00455                 << abort(FatalError);
00456         }
00457     }
00458 
00459     // Correct the face orientation based on the handedness of the block.
00460     // Do nothing for the right-handed block
00461     if (blockHandedness_ == noPoints)
00462     {
00463         FatalErrorIn
00464         (
00465             "patchFaces(const label direc, const labelList& range) const"
00466         )   << "Unable to determine block handedness as points "
00467             << "have not been read in yet"
00468             << abort(FatalError);
00469     }
00470     else if (blockHandedness_ == left)
00471     {
00472         // turn all faces inside out
00473         forAll (result, faceI)
00474         {
00475             result[faceI] = result[faceI].reverseFace();
00476         }
00477     }
00478 
00479     return result;
00480 }
00481 
00482 
00483 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00484 
00485 } // End namespace Foam
00486 
00487 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines