blockCreate.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-2016 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 "error.H"
27 #include "block.H"
28 
29 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
30 
31 void Foam::block::createPoints() const
32 {
33  // Set local variables for mesh specification
34  const label ni = meshDensity().x();
35  const label nj = meshDensity().y();
36  const label nk = meshDensity().z();
37 
38  const point& p000 = blockPoint(0);
39  const point& p100 = blockPoint(1);
40  const point& p110 = blockPoint(2);
41  const point& p010 = blockPoint(3);
42 
43  const point& p001 = blockPoint(4);
44  const point& p101 = blockPoint(5);
45  const point& p111 = blockPoint(6);
46  const point& p011 = blockPoint(7);
47 
48 
49  // List of edge point and weighting factors
50  const List<List<point>>& p = blockEdgePoints();
51  const scalarListList& w = blockEdgeWeights();
52 
53  //
54  // Generate vertices
55  //
56  vertices_.clear();
57  vertices_.setSize(nPoints());
58 
59  for (label k = 0; k <= nk; k++)
60  {
61  for (label j = 0; j <= nj; j++)
62  {
63  for (label i = 0; i <= ni; i++)
64  {
65  const label vertexNo = vtxLabel(i, j, k);
66 
67  // Points on edges
68  vector edgex1 = p000 + (p100 - p000)*w[0][i];
69  vector edgex2 = p010 + (p110 - p010)*w[1][i];
70  vector edgex3 = p011 + (p111 - p011)*w[2][i];
71  vector edgex4 = p001 + (p101 - p001)*w[3][i];
72 
73  vector edgey1 = p000 + (p010 - p000)*w[4][j];
74  vector edgey2 = p100 + (p110 - p100)*w[5][j];
75  vector edgey3 = p101 + (p111 - p101)*w[6][j];
76  vector edgey4 = p001 + (p011 - p001)*w[7][j];
77 
78  vector edgez1 = p000 + (p001 - p000)*w[8][k];
79  vector edgez2 = p100 + (p101 - p100)*w[9][k];
80  vector edgez3 = p110 + (p111 - p110)*w[10][k];
81  vector edgez4 = p010 + (p011 - p010)*w[11][k];
82 
83 
84  // Calculate the importance factors for all edges
85 
86  // x-direction
87  scalar impx1 =
88  (
89  (1.0 - w[0][i])*(1.0 - w[4][j])*(1.0 - w[8][k])
90  + w[0][i]*(1.0 - w[5][j])*(1.0 - w[9][k])
91  );
92 
93  scalar impx2 =
94  (
95  (1.0 - w[1][i])*w[4][j]*(1.0 - w[11][k])
96  + w[1][i]*w[5][j]*(1.0 - w[10][k])
97  );
98 
99  scalar impx3 =
100  (
101  (1.0 - w[2][i])*w[7][j]*w[11][k]
102  + w[2][i]*w[6][j]*w[10][k]
103  );
104 
105  scalar impx4 =
106  (
107  (1.0 - w[3][i])*(1.0 - w[7][j])*w[8][k]
108  + w[3][i]*(1.0 - w[6][j])*w[9][k]
109  );
110 
111  scalar magImpx = impx1 + impx2 + impx3 + impx4;
112  impx1 /= magImpx;
113  impx2 /= magImpx;
114  impx3 /= magImpx;
115  impx4 /= magImpx;
116 
117 
118  // y-direction
119  scalar impy1 =
120  (
121  (1.0 - w[4][j])*(1.0 - w[0][i])*(1.0 - w[8][k])
122  + w[4][j]*(1.0 - w[1][i])*(1.0 - w[11][k])
123  );
124 
125  scalar impy2 =
126  (
127  (1.0 - w[5][j])*w[0][i]*(1.0 - w[9][k])
128  + w[5][j]*w[1][i]*(1.0 - w[10][k])
129  );
130 
131  scalar impy3 =
132  (
133  (1.0 - w[6][j])*w[3][i]*w[9][k]
134  + w[6][j]*w[2][i]*w[10][k]
135  );
136 
137  scalar impy4 =
138  (
139  (1.0 - w[7][j])*(1.0 - w[3][i])*w[8][k]
140  + w[7][j]*(1.0 - w[2][i])*w[11][k]
141  );
142 
143  scalar magImpy = impy1 + impy2 + impy3 + impy4;
144  impy1 /= magImpy;
145  impy2 /= magImpy;
146  impy3 /= magImpy;
147  impy4 /= magImpy;
148 
149 
150  // z-direction
151  scalar impz1 =
152  (
153  (1.0 - w[8][k])*(1.0 - w[0][i])*(1.0 - w[4][j])
154  + w[8][k]*(1.0 - w[3][i])*(1.0 - w[7][j])
155  );
156 
157  scalar impz2 =
158  (
159  (1.0 - w[9][k])*w[0][i]*(1.0 - w[5][j])
160  + w[9][k]*w[3][i]*(1.0 - w[6][j])
161  );
162 
163  scalar impz3 =
164  (
165  (1.0 - w[10][k])*w[1][i]*w[5][j]
166  + w[10][k]*w[2][i]*w[6][j]
167  );
168 
169  scalar impz4 =
170  (
171  (1.0 - w[11][k])*(1.0 - w[1][i])*w[4][j]
172  + w[11][k]*(1.0 - w[2][i])*w[7][j]
173  );
174 
175  scalar magImpz = impz1 + impz2 + impz3 + impz4;
176  impz1 /= magImpz;
177  impz2 /= magImpz;
178  impz3 /= magImpz;
179  impz4 /= magImpz;
180 
181 
182  // Calculate the correction vectors
183  vector corx1 = impx1*(p[0][i] - edgex1);
184  vector corx2 = impx2*(p[1][i] - edgex2);
185  vector corx3 = impx3*(p[2][i] - edgex3);
186  vector corx4 = impx4*(p[3][i] - edgex4);
187 
188  vector cory1 = impy1*(p[4][j] - edgey1);
189  vector cory2 = impy2*(p[5][j] - edgey2);
190  vector cory3 = impy3*(p[6][j] - edgey3);
191  vector cory4 = impy4*(p[7][j] - edgey4);
192 
193  vector corz1 = impz1*(p[8][k] - edgez1);
194  vector corz2 = impz2*(p[9][k] - edgez2);
195  vector corz3 = impz3*(p[10][k] - edgez3);
196  vector corz4 = impz4*(p[11][k] - edgez4);
197 
198 
199  // Multiply by the importance factor
200 
201  // x-direction
202  edgex1 *= impx1;
203  edgex2 *= impx2;
204  edgex3 *= impx3;
205  edgex4 *= impx4;
206 
207  // y-direction
208  edgey1 *= impy1;
209  edgey2 *= impy2;
210  edgey3 *= impy3;
211  edgey4 *= impy4;
212 
213  // z-direction
214  edgez1 *= impz1;
215  edgez2 *= impz2;
216  edgez3 *= impz3;
217  edgez4 *= impz4;
218 
219 
220  // Add the contributions
221  vertices_[vertexNo] =
222  (
223  edgex1 + edgex2 + edgex3 + edgex4
224  + edgey1 + edgey2 + edgey3 + edgey4
225  + edgez1 + edgez2 + edgez3 + edgez4
226  ) / 3.0;
227 
228  vertices_[vertexNo] +=
229  (
230  corx1 + corx2 + corx3 + corx4
231  + cory1 + cory2 + cory3 + cory4
232  + corz1 + corz2 + corz3 + corz4
233  );
234  }
235  }
236  }
237 }
238 
239 
240 void Foam::block::createCells() const
241 {
242  const label ni = meshDensity().x();
243  const label nj = meshDensity().y();
244  const label nk = meshDensity().z();
245 
246  //
247  // Generate cells
248  //
249  cells_.clear();
250  cells_.setSize(nCells());
251 
252  label cellNo = 0;
253 
254  for (label k = 0; k < nk; k++)
255  {
256  for (label j = 0; j < nj; j++)
257  {
258  for (label i = 0; i < ni; i++)
259  {
260  cells_[cellNo].setSize(8);
261 
262  cells_[cellNo][0] = vtxLabel(i, j, k);
263  cells_[cellNo][1] = vtxLabel(i+1, j, k);
264  cells_[cellNo][2] = vtxLabel(i+1, j+1, k);
265  cells_[cellNo][3] = vtxLabel(i, j+1, k);
266  cells_[cellNo][4] = vtxLabel(i, j, k+1);
267  cells_[cellNo][5] = vtxLabel(i+1, j, k+1);
268  cells_[cellNo][6] = vtxLabel(i+1, j+1, k+1);
269  cells_[cellNo][7] = vtxLabel(i, j+1, k+1);
270  cellNo++;
271  }
272  }
273  }
274 }
275 
276 
277 void Foam::block::createBoundary() const
278 {
279  const label ni = meshDensity().x();
280  const label nj = meshDensity().y();
281  const label nk = meshDensity().z();
282 
283  //
284  // Generate boundaries on each side of the hex
285  //
286  boundaryPatches_.clear();
287  boundaryPatches_.setSize(6);
288 
289 
290  // x-direction
291 
292  label wallLabel = 0;
293  label wallCellLabel = 0;
294 
295  // x-min
296  boundaryPatches_[wallLabel].setSize(nj*nk);
297  for (label k = 0; k < nk; k++)
298  {
299  for (label j = 0; j < nj; j++)
300  {
301  boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
302 
303  // Set the points
304  boundaryPatches_[wallLabel][wallCellLabel][0] =
305  vtxLabel(0, j, k);
306  boundaryPatches_[wallLabel][wallCellLabel][1] =
307  vtxLabel(0, j, k + 1);
308  boundaryPatches_[wallLabel][wallCellLabel][2] =
309  vtxLabel(0, j + 1, k + 1);
310  boundaryPatches_[wallLabel][wallCellLabel][3] =
311  vtxLabel(0, j + 1, k);
312 
313  // Update the counter
314  wallCellLabel++;
315  }
316  }
317 
318  // x-max
319  wallLabel++;
320  wallCellLabel = 0;
321 
322  boundaryPatches_[wallLabel].setSize(nj*nk);
323 
324  for (label k = 0; k < nk; k++)
325  {
326  for (label j = 0; j < nj; j++)
327  {
328  boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
329 
330  // Set the points
331  boundaryPatches_[wallLabel][wallCellLabel][0] =
332  vtxLabel(ni, j, k);
333  boundaryPatches_[wallLabel][wallCellLabel][1] =
334  vtxLabel(ni, j+1, k);
335  boundaryPatches_[wallLabel][wallCellLabel][2] =
336  vtxLabel(ni, j+1, k+1);
337  boundaryPatches_[wallLabel][wallCellLabel][3] =
338  vtxLabel(ni, j, k+1);
339 
340  // Update the counter
341  wallCellLabel++;
342  }
343  }
344 
345  // y-direction
346 
347  // y-min
348  wallLabel++;
349  wallCellLabel = 0;
350 
351  boundaryPatches_[wallLabel].setSize(ni*nk);
352  for (label i = 0; i < ni; i++)
353  {
354  for (label k = 0; k < nk; k++)
355  {
356  boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
357 
358  // Set the points
359  boundaryPatches_[wallLabel][wallCellLabel][0] =
360  vtxLabel(i, 0, k);
361  boundaryPatches_[wallLabel][wallCellLabel][1] =
362  vtxLabel(i + 1, 0, k);
363  boundaryPatches_[wallLabel][wallCellLabel][2] =
364  vtxLabel(i + 1, 0, k + 1);
365  boundaryPatches_[wallLabel][wallCellLabel][3] =
366  vtxLabel(i, 0, k + 1);
367 
368  // Update the counter
369  wallCellLabel++;
370  }
371  }
372 
373  // y-max
374  wallLabel++;
375  wallCellLabel = 0;
376 
377  boundaryPatches_[wallLabel].setSize(ni*nk);
378 
379  for (label i = 0; i < ni; i++)
380  {
381  for (label k = 0; k < nk; k++)
382  {
383  boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
384 
385  // Set the points
386  boundaryPatches_[wallLabel][wallCellLabel][0] =
387  vtxLabel(i, nj, k);
388  boundaryPatches_[wallLabel][wallCellLabel][1] =
389  vtxLabel(i, nj, k + 1);
390  boundaryPatches_[wallLabel][wallCellLabel][2] =
391  vtxLabel(i + 1, nj, k + 1);
392  boundaryPatches_[wallLabel][wallCellLabel][3] =
393  vtxLabel(i + 1, nj, k);
394 
395  // Update the counter
396  wallCellLabel++;
397  }
398  }
399 
400  // z-direction
401 
402  // z-min
403  wallLabel++;
404  wallCellLabel = 0;
405 
406  boundaryPatches_[wallLabel].setSize(ni*nj);
407 
408  for (label i = 0; i < ni; i++)
409  {
410  for (label j = 0; j < nj; j++)
411  {
412  boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
413 
414  // Set the points
415  boundaryPatches_[wallLabel][wallCellLabel][0] =
416  vtxLabel(i, j, 0);
417  boundaryPatches_[wallLabel][wallCellLabel][1] =
418  vtxLabel(i, j + 1, 0);
419  boundaryPatches_[wallLabel][wallCellLabel][2] =
420  vtxLabel(i + 1, j + 1, 0);
421  boundaryPatches_[wallLabel][wallCellLabel][3] =
422  vtxLabel(i + 1, j, 0);
423 
424  // Update the counter
425  wallCellLabel++;
426  }
427  }
428 
429  // z-max
430  wallLabel++;
431  wallCellLabel = 0;
432 
433  boundaryPatches_[wallLabel].setSize(ni*nj);
434 
435  for (label i = 0; i < ni; i++)
436  {
437  for (label j = 0; j < nj; j++)
438  {
439  boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
440 
441  // Set the points
442  boundaryPatches_[wallLabel][wallCellLabel][0] =
443  vtxLabel(i, j, nk);
444  boundaryPatches_[wallLabel][wallCellLabel][1] =
445  vtxLabel(i + 1, j, nk);
446  boundaryPatches_[wallLabel][wallCellLabel][2] =
447  vtxLabel(i + 1, j + 1, nk);
448  boundaryPatches_[wallLabel][wallCellLabel][3] =
449  vtxLabel(i, j + 1, nk);
450 
451  // Update the counter
452  wallCellLabel++;
453  }
454  }
455 }
456 
457 
458 void Foam::block::clearGeom()
459 {
460  vertices_.clear();
461  cells_.clear();
462  boundaryPatches_.clear();
463 }
464 
465 
466 // ************************************************************************* //
Foam::Vector::z
const Cmpt & z() const
Definition: VectorI.H:87
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::blockDescriptor::blockEdgeWeights
const scalarListList & blockEdgeWeights() const
Return the weightings along each edge.
Definition: blockDescriptor.C:196
Foam::Vector::x
const Cmpt & x() const
Definition: VectorI.H:75
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
k
label k
Boltzmann constant.
Definition: LISASMDCalcMethod2.H:41
Foam::blockDescriptor::meshDensity
const Vector< label > & meshDensity() const
Return the mesh density (number of cells) in the i,j,k directions.
Definition: blockDescriptor.C:202
Foam::blockDescriptor::blockPoint
const point & blockPoint(const label i) const
Return block point at local label i.
Definition: blockDescriptor.C:236
Foam::block::clearGeom
void clearGeom()
Clear geometry (internal points, cells, boundaryPatches)
Definition: blockCreate.C:458
Foam::blockDescriptor::nPoints
label nPoints() const
Return the number of points.
Definition: blockDescriptor.C:214
Foam::List::clear
void clear()
Clear the list, i.e. set size to zero.
Definition: List.C:356
Foam::scalarListList
List< scalarList > scalarListList
Definition: scalarList.H:51
Foam::Vector::y
const Cmpt & y() const
Definition: VectorI.H:81
Foam::blockDescriptor::nCells
label nCells() const
Return the number of cells.
Definition: blockDescriptor.C:225
Foam::block::vtxLabel
label vtxLabel(label i, label j, label k) const
Vertex label offset for a particular i,j,k position.
Definition: blockI.H:28
Foam::List::setSize
void setSize(const label)
Reset size of List.
Definition: List.C:295
Foam::point
vector point
Point is a vector.
Definition: point.H:41
Foam::blockDescriptor::blockEdgePoints
const List< List< point > > & blockEdgePoints() const
Return the block points along each edge.
Definition: blockDescriptor.C:190
Generated by   doxygen 1.8.11