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 -------------------------------------------------------------------------------
10 License
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 
32 namespace Foam
33 {
34 
35 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
36 
37 label 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)
46 void 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  {
91  WarningInFunction
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
102 hexBlock::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 
114 void 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 
180 labelListList 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  {
238  FatalErrorInFunction
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
253 faceList hexBlock::patchFaces(const label direc, const labelList& range) const
254 {
255  if (range.size() != 6)
256  {
257  FatalErrorInFunction
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  FatalErrorInFunction
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  {
461  FatalErrorInFunction
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 // ************************************************************************* //
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:67
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
nPoints
label nPoints
Definition: gmvOutputHeader.H:2
Foam::constant::physicoChemical::b
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::List::setSize
void setSize(const label n)
Alias for resize()
Definition: List.H:222
Foam::FatalError
error FatalError
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::labelListList
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
range
scalar range
Definition: LISASMDCalcMethod1.H:12
Foam::faceList
List< face > faceList
A List of faces.
Definition: faceListFwd.H:47
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
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.
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
y
scalar y
Definition: LISASMDCalcMethod1.H:14