FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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  WarningIn("hexBlock::hexBlock::setHandedness()")
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 
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 
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  FatalErrorIn("hexBlock::cellShapes()")
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  {
256  (
257  "patchFaces(const label direc, const labelList& range) const"
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  (
453  "patchFaces(const label direc, const labelList& range) const"
454  ) << "direction out of range (1 to 6): " << direc
455  << abort(FatalError);
456  }
457  }
458 
459  // Correct the face orientation based on the handedness of the block.
460  // Do nothing for the right-handed block
461  if (blockHandedness_ == noPoints)
462  {
464  (
465  "patchFaces(const label direc, const labelList& range) const"
466  ) << "Unable to determine block handedness as points "
467  << "have not been read in yet"
468  << abort(FatalError);
469  }
470  else if (blockHandedness_ == left)
471  {
472  // turn all faces inside out
473  forAll (result, faceI)
474  {
475  result[faceI] = result[faceI].reverseFace();
476  }
477  }
478 
479  return result;
480 }
481 
482 
483 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
484 
485 } // End namespace Foam
486 
487 // ************************ vim: set sw=4 sts=4 et: ************************ //