1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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.
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
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
25 \*---------------------------------------------------------------------------*/
29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
34 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
36 label hexBlock::vtxLabel(label a, label b, label c) const
38 return (a + b*(xDim_ + 1) + c*(xDim_ + 1)*(yDim_ + 1));
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
45 void hexBlock::setHandedness()
47 const pointField& p = points_;
49 for (label k = 0; k <= zDim_ - 1; k++)
51 for (label j = 0; j <= yDim_ - 1; j++)
53 for (label i = 0; i <= xDim_ - 1; i++)
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)];
59 if (mag(x) > SMALL && mag(y) > SMALL && mag(z) > SMALL)
61 Info<< "Looking at cell "
62 << i << ' ' << j << ' ' << k
63 << " to determine orientation."
66 if (((x ^ y) & z) > 0)
68 Info<< "Right-handed block." << endl;
69 blockHandedness_ = right;
73 Info << "Left-handed block." << endl;
74 blockHandedness_ = left;
80 Info<< "Cannot determine orientation of cell "
81 << i << ' ' << j << ' ' << k
82 << " since has base vectors " << x << y << z << endl;
88 if (blockHandedness_ == noPoints)
90 WarningIn("hexBlock::hexBlock::setHandedness()")
91 << "Cannot determine orientation of block."
92 << " Continuing as if right handed." << endl;
93 blockHandedness_ = right;
98 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
100 // Construct from components
101 hexBlock::hexBlock(const label nx, const label ny, const label nz)
106 blockHandedness_(noPoints),
107 points_((xDim_ + 1)*(yDim_ + 1)*(zDim_ + 1))
111 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
113 void hexBlock::readPoints
115 const bool readBlank,
116 const scalar twoDThicknes,
122 label nPoints = points_.size();
124 if (twoDThicknes > 0)
129 Info<< "Reading " << nPoints << " x coordinates..." << endl;
130 for (label i=0; i < nPoints; i++)
132 is >> points_[i].x();
135 Info<< "Reading " << nPoints << " y coordinates..." << endl;
136 for (label i=0; i < nPoints; i++)
138 is >> points_[i].y();
141 if (twoDThicknes > 0)
143 Info<< "Extruding " << nPoints << " points in z direction..." << endl;
145 for (label i=0; i < nPoints; i++)
147 points_[i+nPoints] = points_[i];
149 for (label i=0; i < nPoints; i++)
152 points_[i+nPoints].z() = twoDThicknes;
157 Info<< "Reading " << nPoints << " z coordinates..." << endl;
158 for (label i=0; i < nPoints; i++)
160 is >> points_[i].z();
167 Info<< "Reading " << nPoints << " blanks..." << endl;
168 for (label i=0; i < nPoints; i++)
174 // Set left- or righthandedness of block by looking at a cell.
179 labelListList hexBlock::blockCells() const
181 labelListList result(xDim_*yDim_*zDim_);
185 if (blockHandedness_ == right)
187 for (label k = 0; k <= zDim_ - 1; k++)
189 for (label j = 0; j <= yDim_ - 1; j++)
191 for (label i = 0; i <= xDim_ - 1; i++)
193 labelList& hexLabels = result[cellNo];
194 hexLabels.setSize(8);
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);
210 else if (blockHandedness_ == left)
212 for (label k = 0; k <= zDim_ - 1; k++)
214 for (label j = 0; j <= yDim_ - 1; j++)
216 for (label i = 0; i <= xDim_ - 1; i++)
218 labelList& hexLabels = result[cellNo];
219 hexLabels.setSize(8);
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);
237 FatalErrorIn("hexBlock::cellShapes()")
238 << "Unable to determine block handedness as points "
239 << "have not been read in yet"
240 << abort(FatalError);
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
254 if (range.size() != 6)
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);
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];
281 (yMaxRange - yMinRange + 1)*(zMaxRange - zMinRange + 1)
285 for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
287 for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
289 result[p].setSize(4);
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);
310 (xMaxRange - xMinRange + 1)*(zMaxRange - zMinRange + 1)
314 for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
316 for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
318 result[p].setSize(4);
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);
339 (xMaxRange - xMinRange + 1)*(yMaxRange - yMinRange + 1)
343 for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
345 for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
347 result[p].setSize(4);
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_);
368 (yMaxRange - yMinRange + 1)*(zMaxRange - zMinRange + 1)
372 for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
374 for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
376 result[p].setSize(4);
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);
397 (xMaxRange - xMinRange + 1)*(zMaxRange - zMinRange + 1)
401 for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
403 for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
405 result[p].setSize(4);
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);
426 (xMaxRange - xMinRange + 1)*(yMaxRange - yMinRange + 1)
430 for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
432 for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
434 result[p].setSize(4);
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);
454 "patchFaces(const label direc, const labelList& range) const"
455 ) << "direction out of range (1 to 6): " << direc
456 << abort(FatalError);
460 // Correct the face orientation based on the handedness of the block.
461 // Do nothing for the right-handed block
462 if (blockHandedness_ == noPoints)
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);
471 else if (blockHandedness_ == left)
473 // turn all faces inside out
474 forAll (result, faceI)
476 result[faceI] = result[faceI].reverseFace();
484 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
486 } // End namespace Foam
488 // ************************ vim: set sw=4 sts=4 et: ************************ //