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 | Copyright (C) 2011-2015 OpenFOAM Foundation
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
13  the Free Software Foundation, either version 3 of the License, or
14  (at your 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, 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  WarningInFunction
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  FatalErrorInFunction
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  FatalErrorInFunction
256  << "Invalid size of the range array: " << range.size()
257  << ". Should be 6 (xMin, xMax, yMin, yMax, zMin, zMax"
258  << abort(FatalError);
259  }
260 
261  label xMinRange = range[0];
262  label xMaxRange = range[1];
263  label yMinRange = range[2];
264  label yMaxRange = range[3];
265  label zMinRange = range[4];
266  label zMaxRange = range[5];
267 
268  faceList result(0);
269 
270  switch (direc)
271  {
272  case 1:
273  {
274  // high i = xmax
275 
276  result.setSize
277  (
278  (yMaxRange - yMinRange + 1)*(zMaxRange - zMinRange + 1)
279  );
280 
281  label p = 0;
282  for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
283  {
284  for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
285  {
286  result[p].setSize(4);
287 
288  // set the points
289  result[p][0] = vtxLabel(xDim_, j, k);
290  result[p][1] = vtxLabel(xDim_, j+1, k);
291  result[p][2] = vtxLabel(xDim_, j+1, k+1);
292  result[p][3] = vtxLabel(xDim_, j, k+1);
293 
294  p++;
295  }
296  }
297 
298  result.setSize(p);
299  break;
300  }
301 
302  case 2:
303  {
304  // high j = ymax
305  result.setSize
306  (
307  (xMaxRange - xMinRange + 1)*(zMaxRange - zMinRange + 1)
308  );
309 
310  label p = 0;
311  for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
312  {
313  for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
314  {
315  result[p].setSize(4);
316 
317  // set the points
318  result[p][0] = vtxLabel(i, yDim_, k);
319  result[p][1] = vtxLabel(i, yDim_, k + 1);
320  result[p][2] = vtxLabel(i + 1, yDim_, k + 1);
321  result[p][3] = vtxLabel(i + 1, yDim_, k);
322 
323  p++;
324  }
325  }
326 
327  result.setSize(p);
328  break;
329  }
330 
331  case 3:
332  {
333  // high k = zmax
334  result.setSize
335  (
336  (xMaxRange - xMinRange + 1)*(yMaxRange - yMinRange + 1)
337  );
338 
339  label p = 0;
340  for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
341  {
342  for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
343  {
344  result[p].setSize(4);
345 
346  // set the points
347  result[p][0] = vtxLabel(i, j, zDim_);
348  result[p][1] = vtxLabel(i + 1, j, zDim_);
349  result[p][2] = vtxLabel(i + 1, j + 1, zDim_);
350  result[p][3] = vtxLabel(i, j + 1, zDim_);
351 
352  p++;
353  }
354  }
355 
356  result.setSize(p);
357  break;
358  }
359 
360  case 4:
361  {
362  // low i = xmin
363  result.setSize
364  (
365  (yMaxRange - yMinRange + 1)*(zMaxRange - zMinRange + 1)
366  );
367 
368  label p = 0;
369  for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
370  {
371  for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
372  {
373  result[p].setSize(4);
374 
375  // set the points
376  result[p][0] = vtxLabel(0, j, k);
377  result[p][1] = vtxLabel(0, j, k + 1);
378  result[p][2] = vtxLabel(0, j + 1, k + 1);
379  result[p][3] = vtxLabel(0, j + 1, k);
380 
381  p++;
382  }
383  }
384 
385  result.setSize(p);
386  break;
387  }
388 
389  case 5:
390  {
391  // low j = ymin
392  result.setSize
393  (
394  (xMaxRange - xMinRange + 1)*(zMaxRange - zMinRange + 1)
395  );
396 
397  label p = 0;
398  for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
399  {
400  for (label k = zMinRange - 1; k <= zMaxRange - 1; k++)
401  {
402  result[p].setSize(4);
403 
404  // set the points
405  result[p][0] = vtxLabel(i, 0, k);
406  result[p][1] = vtxLabel(i + 1, 0, k);
407  result[p][2] = vtxLabel(i + 1, 0, k + 1);
408  result[p][3] = vtxLabel(i, 0, k + 1);
409 
410  p++;
411  }
412  }
413 
414  result.setSize(p);
415  break;
416  }
417 
418  case 6:
419  {
420  // low k = zmin
421  result.setSize
422  (
423  (xMaxRange - xMinRange + 1)*(yMaxRange - yMinRange + 1)
424  );
425 
426  label p = 0;
427  for (label i = xMinRange - 1; i <= xMaxRange - 1; i++)
428  {
429  for (label j = yMinRange - 1; j <= yMaxRange - 1; j++)
430  {
431  result[p].setSize(4);
432 
433  // set the points
434  result[p][0] = vtxLabel(i, j, 0);
435  result[p][1] = vtxLabel(i, j + 1, 0);
436  result[p][2] = vtxLabel(i + 1, j + 1, 0);
437  result[p][3] = vtxLabel(i + 1, j, 0);
438 
439  p++;
440  }
441  }
442 
443  result.setSize(p);
444  break;
445  }
446 
447  default:
448  {
449  FatalErrorInFunction
450  << "direction out of range (1 to 6): " << direc
451  << abort(FatalError);
452  }
453  }
454 
455  // Correct the face orientation based on the handedness of the block.
456  // Do nothing for the right-handed block
457  if (blockHandedness_ == noPoints)
458  {
459  FatalErrorInFunction
460  << "Unable to determine block handedness as points "
461  << "have not been read in yet"
462  << abort(FatalError);
463  }
464  else if (blockHandedness_ == left)
465  {
466  // turn all faces inside out
467  forAll(result, faceI)
468  {
469  result[faceI].flip();
470  }
471  }
472 
473  return result;
474 }
475 
476 
477 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
478 
479 } // End namespace Foam
480 
481 // ************************************************************************* //
Foam::hexBlock::zDim_
label zDim_
Definition: hexBlock.H:65
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Foam::hexBlock::right
@ right
Definition: hexBlock.H:58
p
p
Definition: pEqn.H:62
Foam::labelList
List< label > labelList
A List of labels.
Definition: labelList.H:56
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:406
Foam::hexBlock::blockCells
labelListList blockCells() const
Return block cells.
Foam::hexBlock::setHandedness
void setHandedness()
Calculate handedness of block.
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Foam::hexBlock::hexBlock
hexBlock(const hexBlock &)
Disallow default bitwise copy construct.
Foam::mag
dimensioned< scalar > mag(const dimensioned< Type > &)
Foam::hexBlock::patchFaces
faceList patchFaces(label direc, const labelList &range) const
Return block patch faces given direction and range limits.
nPoints
label nPoints
Definition: gmvOutputHeader.H:2
Foam::hexBlock::xDim_
label xDim_
Number of point in each direction.
Definition: hexBlock.H:63
Foam::hexBlock::points_
pointField points_
List of points.
Definition: hexBlock.H:71
Foam::label
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Foam::hexBlock::vtxLabel
label vtxLabel(label i, label j, label k) const
Vertex addressing inside the block.
Foam::constant::physicoChemical::b
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:28
Foam::Info
messageStream Info
Foam::hexBlock::left
@ left
Definition: hexBlock.H:59
Foam::FatalError
error FatalError
Foam::hexBlock::noPoints
@ noPoints
Definition: hexBlock.H:57
Foam::hexBlock::yDim_
label yDim_
Definition: hexBlock.H:64
Foam
Namespace for OpenFOAM.
Definition: combustionModel.C:30
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Foam::List::setSize
void setSize(const label)
Reset size of List.
Foam::labelListList
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:318
range
scalar range
Definition: LISASMDCalcMethod1.H:12
Foam::faceList
List< face > faceList
Definition: faceListFwd.H:43
k
label k
Boltzmann constant.
Definition: LISASMDCalcMethod2.H:41
x
x
Definition: LISASMDCalcMethod2.H:52
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
Foam::hexBlock::blockHandedness_
handed blockHandedness_
Handedness of the block.
Definition: hexBlock.H:68
Foam::hexBlock::readPoints
void readPoints(Istream &)
Read block points.
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:259
y
scalar y
Definition: LISASMDCalcMethod1.H:14

Copyright © 2011-2018  OpenFOAM | OPENFOAM® is a registered trademark of OpenCFD Ltd.

Generated by   doxygen 1.8.17