hexBlock.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "hexBlock.H"
29
30// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31
32namespace Foam
33{
34
35// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
36
37label hexBlock::vtxLabel(label a, label b, label c) const
38{
39 return (a + b*(xDim_ + 1) + c*(xDim_ + 1)*(yDim_ + 1));
40}
41
42
43// Calculate the handedness of the block by looking at the orientation
44// of the spanning edges of a hex. Loops until valid cell found (since might
45// be prism)
46void hexBlock::setHandedness()
47{
48 const pointField& p = points_;
49
50 for (label k = 0; k <= zDim_ - 1; k++)
51 {
52 for (label j = 0; j <= yDim_ - 1; j++)
53 {
54 for (label i = 0; i <= xDim_ - 1; i++)
55 {
56 vector x = p[vtxLabel(i+1, j, k)] - p[vtxLabel(i, j, k)];
57 vector y = p[vtxLabel(i, j+1, k)] - p[vtxLabel(i, j, k)];
58 vector z = p[vtxLabel(i, j, k+1)] - p[vtxLabel(i, j, k)];
59
60 if (mag(x) > SMALL && mag(y) > SMALL && mag(z) > SMALL)
61 {
62 Info<< "Looking at cell "
63 << i << ' ' << j << ' ' << k
64 << " to determine orientation."
65 << endl;
66
67 if (((x ^ y) & z) > 0)
68 {
69 Info<< "Right-handed block." << endl;
70 blockHandedness_ = right;
71 }
72 else
73 {
74 Info<< "Left-handed block." << endl;
75 blockHandedness_ = left;
76 }
77 return;
78 }
79 else
80 {
81 Info<< "Cannot determine orientation of cell "
82 << i << ' ' << j << ' ' << k
83 << " since has base vectors " << x << y << z << endl;
84 }
85 }
86 }
87 }
88
89 if (blockHandedness_ == noPoints)
90 {
92 << "Cannot determine orientation of block."
93 << " Continuing as if right handed." << endl;
94 blockHandedness_ = right;
95 }
96}
97
98
99// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
100
101// Construct from components
102hexBlock::hexBlock(const label nx, const label ny, const label nz)
103:
104 xDim_(nx - 1),
105 yDim_(ny - 1),
106 zDim_(nz - 1),
107 blockHandedness_(noPoints),
108 points_((xDim_ + 1)*(yDim_ + 1)*(zDim_ + 1))
109{}
110
111
112// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
113
114void hexBlock::readPoints
115(
116 const bool readBlank,
117 const scalar twoDThicknes,
118 Istream& is
119)
120{
121 scalar iBlank;
122
123 label nPoints = points_.size();
124
125 if (twoDThicknes > 0)
126 {
127 nPoints /= 2;
128 }
129
130 Info<< "Reading " << nPoints << " x coordinates..." << endl;
131 for (label i=0; i < nPoints; i++)
132 {
133 is >> points_[i].x();
134 }
135
136 Info<< "Reading " << nPoints << " y coordinates..." << endl;
137 for (label i=0; i < nPoints; i++)
138 {
139 is >> points_[i].y();
140 }
141
142 if (twoDThicknes > 0)
143 {
144 Info<< "Extruding " << nPoints << " points in z direction..." << endl;
145 // Duplicate points
146 for (label i=0; i < nPoints; i++)
147 {
148 points_[i+nPoints] = points_[i];
149 }
150 for (label i=0; i < nPoints; i++)
151 {
152 points_[i].z() = 0;
153 points_[i+nPoints].z() = twoDThicknes;
154 }
155 }
156 else
157 {
158 Info<< "Reading " << nPoints << " z coordinates..." << endl;
159 for (label i=0; i < nPoints; i++)
160 {
161 is >> points_[i].z();
162 }
163 }
164
165
166 if (readBlank)
167 {
168 Info<< "Reading " << nPoints << " blanks..." << endl;
169 for (label i=0; i < nPoints; i++)
170 {
171 is >> iBlank;
172 }
173 }
174
175 // Set left- or righthandedness of block by looking at a cell.
176 setHandedness();
177}
178
179
180labelListList hexBlock::blockCells() const
181{
182 labelListList result(xDim_*yDim_*zDim_);
183
184 label cellNo = 0;
185
186 if (blockHandedness_ == right)
187 {
188 for (label k = 0; k <= zDim_ - 1; k++)
189 {
190 for (label j = 0; j <= yDim_ - 1; j++)
191 {
192 for (label i = 0; i <= xDim_ - 1; i++)
193 {
194 labelList& hexLabels = result[cellNo];
195 hexLabels.setSize(8);
196
197 hexLabels[0] = vtxLabel(i, j, k);
198 hexLabels[1] = vtxLabel(i+1, j, k);
199 hexLabels[2] = vtxLabel(i+1, j+1, k);
200 hexLabels[3] = vtxLabel(i, j+1, k);
201 hexLabels[4] = vtxLabel(i, j, k+1);
202 hexLabels[5] = vtxLabel(i+1, j, k+1);
203 hexLabels[6] = vtxLabel(i+1, j+1, k+1);
204 hexLabels[7] = vtxLabel(i, j+1, k+1);
205
206 cellNo++;
207 }
208 }
209 }
210 }
211 else if (blockHandedness_ == left)
212 {
213 for (label k = 0; k <= zDim_ - 1; k++)
214 {
215 for (label j = 0; j <= yDim_ - 1; j++)
216 {
217 for (label i = 0; i <= xDim_ - 1; i++)
218 {
219 labelList& hexLabels = result[cellNo];
220 hexLabels.setSize(8);
221
222 hexLabels[0] = vtxLabel(i, j, k+1);
223 hexLabels[1] = vtxLabel(i+1, j, k+1);
224 hexLabels[2] = vtxLabel(i+1, j+1, k+1);
225 hexLabels[3] = vtxLabel(i, j+1, k+1);
226 hexLabels[4] = vtxLabel(i, j, k);
227 hexLabels[5] = vtxLabel(i+1, j, k);
228 hexLabels[6] = vtxLabel(i+1, j+1, k);
229 hexLabels[7] = vtxLabel(i, j+1, k);
230
231 cellNo++;
232 }
233 }
234 }
235 }
236 else
237 {
239 << "Unable to determine block handedness as points "
240 << "have not been read in yet"
241 << abort(FatalError);
242 }
243
244 return result;
245}
246
247
248// Return block patch faces given direction and range limits
249// From the cfx manual: direction
250// 0 = solid (3-D patch),
251// 1 = high i, 2 = high j, 3 = high k
252// 4 = low i, 5 = low j, 6 = low k
253faceList hexBlock::patchFaces(const label direc, const labelList& range) const
254{
255 if (range.size() != 6)
256 {
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 {
452 << "direction out of range (1 to 6): " << direc
453 << abort(FatalError);
454 }
455 }
456
457 // Correct the face orientation based on the handedness of the block.
458 // Do nothing for the right-handed block
459 if (blockHandedness_ == noPoints)
460 {
462 << "Unable to determine block handedness as points "
463 << "have not been read in yet"
464 << abort(FatalError);
465 }
466 else if (blockHandedness_ == left)
467 {
468 // turn all faces inside out
469 forAll(result, facei)
470 {
471 result[facei].flip();
472 }
473 }
474
475 return result;
476}
477
478
479// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
480
481} // End namespace Foam
482
483// ************************************************************************* //
scalar range
scalar y
label k
void setSize(const label n)
Alias for resize()
Definition: List.H:218
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
volScalarField & p
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
label nPoints
#define WarningInFunction
Report a warning using Foam::Warning.
Namespace for OpenFOAM.
List< label > labelList
A List of labels.
Definition: List.H:66
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
errorManip< error > abort(error &err)
Definition: errorManip.H:144
List< face > faceList
A List of faces.
Definition: faceListFwd.H:47
volScalarField & b
Definition: createFields.H:27
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333