blockCreate.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 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 "block.H"
27 
28 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
29 
30 #define w0 w[0][i]
31 #define w1 w[1][i]
32 #define w2 w[2][i]
33 #define w3 w[3][i]
34 
35 #define w4 w[4][j]
36 #define w5 w[5][j]
37 #define w6 w[6][j]
38 #define w7 w[7][j]
39 
40 #define w8 w[8][k]
41 #define w9 w[9][k]
42 #define w10 w[10][k]
43 #define w11 w[11][k]
44 
45 void Foam::block::createPoints()
46 {
47  // Set local variables for mesh specification
48  const label ni = density().x();
49  const label nj = density().y();
50  const label nk = density().z();
51 
52  const point& p000 = blockPoint(0);
53  const point& p100 = blockPoint(1);
54  const point& p110 = blockPoint(2);
55  const point& p010 = blockPoint(3);
56 
57  const point& p001 = blockPoint(4);
58  const point& p101 = blockPoint(5);
59  const point& p111 = blockPoint(6);
60  const point& p011 = blockPoint(7);
61 
62  // List of edge point and weighting factors
63  pointField p[12];
64  scalarList w[12];
65  label nCurvedEdges = edgesPointsWeights(p, w);
66 
67  points_.setSize(nPoints());
68 
69  points_[pointLabel(0, 0, 0)] = p000;
70  points_[pointLabel(ni, 0, 0)] = p100;
71  points_[pointLabel(ni, nj, 0)] = p110;
72  points_[pointLabel(0, nj, 0)] = p010;
73  points_[pointLabel(0, 0, nk)] = p001;
74  points_[pointLabel(ni, 0, nk)] = p101;
75  points_[pointLabel(ni, nj, nk)] = p111;
76  points_[pointLabel(0, nj, nk)] = p011;
77 
78  for (label k=0; k<=nk; k++)
79  {
80  for (label j=0; j<=nj; j++)
81  {
82  for (label i=0; i<=ni; i++)
83  {
84  // Skip block vertices
85  if (vertex(i, j, k)) continue;
86 
87  const label vijk = pointLabel(i, j, k);
88 
89  // Calculate the weighting factors for all edges
90 
91  // x-direction
92  scalar wx1 = (1 - w0)*(1 - w4)*(1 - w8) + w0*(1 - w5)*(1 - w9);
93  scalar wx2 = (1 - w1)*w4*(1 - w11) + w1*w5*(1 - w10);
94  scalar wx3 = (1 - w2)*w7*w11 + w2*w6*w10;
95  scalar wx4 = (1 - w3)*(1 - w7)*w8 + w3*(1 - w6)*w9;
96 
97  const scalar sumWx = wx1 + wx2 + wx3 + wx4;
98  wx1 /= sumWx;
99  wx2 /= sumWx;
100  wx3 /= sumWx;
101  wx4 /= sumWx;
102 
103 
104  // y-direction
105  scalar wy1 = (1 - w4)*(1 - w0)*(1 - w8) + w4*(1 - w1)*(1 - w11);
106  scalar wy2 = (1 - w5)*w0*(1 - w9) + w5*w1*(1 - w10);
107  scalar wy3 = (1 - w6)*w3*w9 + w6*w2*w10;
108  scalar wy4 = (1 - w7)*(1 - w3)*w8 + w7*(1 - w2)*w11;
109 
110  const scalar sumWy = wy1 + wy2 + wy3 + wy4;
111  wy1 /= sumWy;
112  wy2 /= sumWy;
113  wy3 /= sumWy;
114  wy4 /= sumWy;
115 
116 
117  // z-direction
118  scalar wz1 = (1 - w8)*(1 - w0)*(1 - w4) + w8*(1 - w3)*(1 - w7);
119  scalar wz2 = (1 - w9)*w0*(1 - w5) + w9*w3*(1 - w6);
120  scalar wz3 = (1 - w10)*w1*w5 + w10*w2*w6;
121  scalar wz4 = (1 - w11)*(1 - w1)*w4 + w11*(1 - w2)*w7;
122 
123  const scalar sumWz = wz1 + wz2 + wz3 + wz4;
124  wz1 /= sumWz;
125  wz2 /= sumWz;
126  wz3 /= sumWz;
127  wz4 /= sumWz;
128 
129 
130  // Points on straight edges
131  const vector edgex1 = p000 + (p100 - p000)*w0;
132  const vector edgex2 = p010 + (p110 - p010)*w1;
133  const vector edgex3 = p011 + (p111 - p011)*w2;
134  const vector edgex4 = p001 + (p101 - p001)*w3;
135 
136  const vector edgey1 = p000 + (p010 - p000)*w4;
137  const vector edgey2 = p100 + (p110 - p100)*w5;
138  const vector edgey3 = p101 + (p111 - p101)*w6;
139  const vector edgey4 = p001 + (p011 - p001)*w7;
140 
141  const vector edgez1 = p000 + (p001 - p000)*w8;
142  const vector edgez2 = p100 + (p101 - p100)*w9;
143  const vector edgez3 = p110 + (p111 - p110)*w10;
144  const vector edgez4 = p010 + (p011 - p010)*w11;
145 
146  // Add the contributions
147  points_[vijk] =
148  (
149  wx1*edgex1 + wx2*edgex2 + wx3*edgex3 + wx4*edgex4
150  + wy1*edgey1 + wy2*edgey2 + wy3*edgey3 + wy4*edgey4
151  + wz1*edgez1 + wz2*edgez2 + wz3*edgez3 + wz4*edgez4
152  )/3;
153 
154 
155  // Apply curved-edge correction if block has curved edges
156  if (nCurvedEdges)
157  {
158  // Calculate the correction vectors
159  const vector corx1 = wx1*(p[0][i] - edgex1);
160  const vector corx2 = wx2*(p[1][i] - edgex2);
161  const vector corx3 = wx3*(p[2][i] - edgex3);
162  const vector corx4 = wx4*(p[3][i] - edgex4);
163 
164  const vector cory1 = wy1*(p[4][j] - edgey1);
165  const vector cory2 = wy2*(p[5][j] - edgey2);
166  const vector cory3 = wy3*(p[6][j] - edgey3);
167  const vector cory4 = wy4*(p[7][j] - edgey4);
168 
169  const vector corz1 = wz1*(p[8][k] - edgez1);
170  const vector corz2 = wz2*(p[9][k] - edgez2);
171  const vector corz3 = wz3*(p[10][k] - edgez3);
172  const vector corz4 = wz4*(p[11][k] - edgez4);
173 
174  points_[vijk] +=
175  (
176  corx1 + corx2 + corx3 + corx4
177  + cory1 + cory2 + cory3 + cory4
178  + corz1 + corz2 + corz3 + corz4
179  );
180  }
181  }
182  }
183  }
184 
185  if (!nCurvedFaces()) return;
186 
187  // Apply curvature correction to face points
188  FixedList<pointField, 6> facePoints(this->facePoints(points_));
189  correctFacePoints(facePoints);
190 
191  // Loop over points and apply the correction from the from the i-faces
192  for (label ii=0; ii<=ni; ii++)
193  {
194  // Process the points on the curved faces last
195  label i = (ii + 1)%(ni + 1);
196 
197  for (label j=0; j<=nj; j++)
198  {
199  for (label k=0; k<=nk; k++)
200  {
201  // Skip non-curved faces and edges
202  if (flatFaceOrEdge(i, j, k)) continue;
203 
204  const label vijk = pointLabel(i, j, k);
205 
206  scalar wf0 =
207  (
208  (1 - w0)*(1 - w4)*(1 - w8)
209  + (1 - w1)*w4*(1 - w11)
210  + (1 - w2)*w7*w11
211  + (1 - w3)*(1 - w7)*w8
212  );
213 
214  scalar wf1 =
215  (
216  w0*(1 - w5)*(1 - w9)
217  + w1*w5*(1 - w10)
218  + w2*w5*w10
219  + w3*(1 - w6)*w9
220  );
221 
222  const scalar sumWf = wf0 + wf1;
223  wf0 /= sumWf;
224  wf1 /= sumWf;
225 
226  points_[vijk] +=
227  (
228  wf0
229  *(
230  facePoints[0][facePointLabel(0, j, k)]
231  - points_[pointLabel(0, j, k)]
232  )
233  + wf1
234  *(
235  facePoints[1][facePointLabel(1, j, k)]
236  - points_[pointLabel(ni, j, k)]
237  )
238  );
239  }
240  }
241  }
242 
243  // Loop over points and apply the correction from the from the j-faces
244  for (label jj=0; jj<=nj; jj++)
245  {
246  // Process the points on the curved faces last
247  label j = (jj + 1)%(nj + 1);
248 
249  for (label i=0; i<=ni; i++)
250  {
251  for (label k=0; k<=nk; k++)
252  {
253  // Skip flat faces and edges
254  if (flatFaceOrEdge(i, j, k)) continue;
255 
256  const label vijk = pointLabel(i, j, k);
257 
258  scalar wf2 =
259  (
260  (1 - w4)*(1 - w1)*(1 - w8)
261  + (1 - w5)*w0*(1 - w9)
262  + (1 - w6)*w3*w9
263  + (1 - w7)*(1 - w3)*w8
264  );
265 
266  scalar wf3 =
267  (
268  w4*(1 - w1)*(1 - w11)
269  + w5*w1*(1 - w10)
270  + w6*w2*w10
271  + w7*(1 - w2)*w11
272  );
273 
274  const scalar sumWf = wf2 + wf3;
275  wf2 /= sumWf;
276  wf3 /= sumWf;
277 
278  points_[vijk] +=
279  (
280  wf2
281  *(
282  facePoints[2][facePointLabel(2, i, k)]
283  - points_[pointLabel(i, 0, k)]
284  )
285  + wf3
286  *(
287  facePoints[3][facePointLabel(3, i, k)]
288  - points_[pointLabel(i, nj, k)]
289  )
290  );
291  }
292  }
293  }
294 
295  // Loop over points and apply the correction from the from the k-faces
296  for (label kk=0; kk<=nk; kk++)
297  {
298  // Process the points on the curved faces last
299  label k = (kk + 1)%(nk + 1);
300 
301  for (label i=0; i<=ni; i++)
302  {
303  for (label j=0; j<=nj; j++)
304  {
305  // Skip flat faces and edges
306  if (flatFaceOrEdge(i, j, k)) continue;
307 
308  const label vijk = pointLabel(i, j, k);
309 
310  scalar wf4 =
311  (
312  (1 - w8)*(1 - w0)*(1 - w4)
313  + (1 - w9)*w0*(1 - w5)
314  + (1 - w10)*w1*w5
315  + (1 - w11)*(1 - w1)*w4
316  );
317 
318  scalar wf5 =
319  (
320  w8*(1 - w3)*(1 - w7)
321  + w9*w3*(1 - w6)
322  + w10*w2*w6
323  + w11*(1 - w2)*w7
324  );
325 
326  const scalar sumWf = wf4 + wf5;
327  wf4 /= sumWf;
328  wf5 /= sumWf;
329 
330  points_[vijk] +=
331  (
332  wf4
333  *(
334  facePoints[4][facePointLabel(4, i, j)]
335  - points_[pointLabel(i, j, 0)]
336  )
337  + wf5
338  *(
339  facePoints[5][facePointLabel(5, i, j)]
340  - points_[pointLabel(i, j, nk)]
341  )
342  );
343  }
344  }
345  }
346 }
347 
348 
350 {
351  const label ni = density().x();
352  const label nj = density().y();
353  const label nk = density().z();
354 
356 
357  label celli = 0;
358 
359  for (label k=0; k<nk; k++)
360  {
361  for (label j=0; j<nj; j++)
362  {
363  for (label i=0; i<ni; i++)
364  {
365  cells[celli][0] = pointLabel(i, j, k);
366  cells[celli][1] = pointLabel(i+1, j, k);
367  cells[celli][2] = pointLabel(i+1, j+1, k);
368  cells[celli][3] = pointLabel(i, j+1, k);
369  cells[celli][4] = pointLabel(i, j, k+1);
370  cells[celli][5] = pointLabel(i+1, j, k+1);
371  cells[celli][6] = pointLabel(i+1, j+1, k+1);
372  cells[celli][7] = pointLabel(i, j+1, k+1);
373 
374  celli++;
375  }
376  }
377  }
378 
379  return cells;
380 }
381 
382 
383 void Foam::block::createBoundary()
384 {
385  const label ni = density().x();
386  const label nj = density().y();
387  const label nk = density().z();
388 
389  label patchi = 0;
390  label facei = 0;
391 
392  // x-direction
393 
394  // x-min
395  boundaryPatches_[patchi].setSize(nj*nk);
396  for (label k=0; k<nk; k++)
397  {
398  for (label j=0; j<nj; j++)
399  {
400  boundaryPatches_[patchi][facei][0] = pointLabel(0, j, k);
401  boundaryPatches_[patchi][facei][1] = pointLabel(0, j, k+1);
402  boundaryPatches_[patchi][facei][2] = pointLabel(0, j+1, k+1);
403  boundaryPatches_[patchi][facei][3] = pointLabel(0, j+1, k);
404 
405  facei++;
406  }
407  }
408 
409  // x-max
410  patchi++;
411  facei = 0;
412 
413  boundaryPatches_[patchi].setSize(nj*nk);
414 
415  for (label k=0; k<nk; k++)
416  {
417  for (label j=0; j<nj; j++)
418  {
419  boundaryPatches_[patchi][facei][0] = pointLabel(ni, j, k);
420  boundaryPatches_[patchi][facei][1] = pointLabel(ni, j+1, k);
421  boundaryPatches_[patchi][facei][2] = pointLabel(ni, j+1, k+1);
422  boundaryPatches_[patchi][facei][3] = pointLabel(ni, j, k+1);
423 
424  facei++;
425  }
426  }
427 
428  // y-direction
429 
430  // y-min
431  patchi++;
432  facei = 0;
433 
434  boundaryPatches_[patchi].setSize(ni*nk);
435  for (label i=0; i<ni; i++)
436  {
437  for (label k=0; k<nk; k++)
438  {
439  boundaryPatches_[patchi][facei][0] = pointLabel(i, 0, k);
440  boundaryPatches_[patchi][facei][1] = pointLabel(i+1, 0, k);
441  boundaryPatches_[patchi][facei][2] = pointLabel(i+1, 0, k+1);
442  boundaryPatches_[patchi][facei][3] = pointLabel(i, 0, k+1);
443 
444  facei++;
445  }
446  }
447 
448  // y-max
449  patchi++;
450  facei = 0;
451 
452  boundaryPatches_[patchi].setSize(ni*nk);
453 
454  for (label i=0; i<ni; i++)
455  {
456  for (label k=0; k<nk; k++)
457  {
458  boundaryPatches_[patchi][facei][0] = pointLabel(i, nj, k);
459  boundaryPatches_[patchi][facei][1] = pointLabel(i, nj, k+1);
460  boundaryPatches_[patchi][facei][2] = pointLabel(i+1, nj, k+1);
461  boundaryPatches_[patchi][facei][3] = pointLabel(i+1, nj, k);
462 
463  facei++;
464  }
465  }
466 
467  // z-direction
468 
469  // z-min
470  patchi++;
471  facei = 0;
472 
473  boundaryPatches_[patchi].setSize(ni*nj);
474 
475  for (label i=0; i<ni; i++)
476  {
477  for (label j=0; j<nj; j++)
478  {
479  boundaryPatches_[patchi][facei][0] = pointLabel(i, j, 0);
480  boundaryPatches_[patchi][facei][1] = pointLabel(i, j+1, 0);
481  boundaryPatches_[patchi][facei][2] = pointLabel(i+1, j+1, 0);
482  boundaryPatches_[patchi][facei][3] = pointLabel(i+1, j, 0);
483 
484  facei++;
485  }
486  }
487 
488  // z-max
489  patchi++;
490  facei = 0;
491 
492  boundaryPatches_[patchi].setSize(ni*nj);
493 
494  for (label i=0; i<ni; i++)
495  {
496  for (label j=0; j<nj; j++)
497  {
498  boundaryPatches_[patchi][facei][0] = pointLabel(i, j, nk);
499  boundaryPatches_[patchi][facei][1] = pointLabel(i+1, j, nk);
500  boundaryPatches_[patchi][facei][2] = pointLabel(i+1, j+1, nk);
501  boundaryPatches_[patchi][facei][3] = pointLabel(i, j+1, nk);
502 
503  facei++;
504  }
505  }
506 }
507 
508 
509 // ************************************************************************* //
label edgesPointsWeights(pointField(&edgePoints)[12], scalarList(&edgeWeights)[12]) const
Calculate the points and weights for all edges.
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
bool vertex(const label i, const label j, const label k) const
Return true if point i,j,k addresses a block vertex.
const Vector< label > & density() const
Return the mesh density (number of cells) in the i,j,k directions.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
label nCells() const
Return the number of cells.
List< FixedList< label, 8 > > cells() const
Return the cells for filling the block.
Definition: blockCreate.C:349
const Cmpt & z() const
Definition: VectorI.H:87
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
#define w2
Definition: blockCreate.C:32
label k
Boltzmann constant.
FixedList< pointField, 6 > facePoints(const pointField &points) const
Return the list of face-points for all of the faces of the block.
#define w4
Definition: blockCreate.C:35
#define w11
Definition: blockCreate.C:43
const Cmpt & y() const
Definition: VectorI.H:81
const point & blockPoint(const label i) const
Return block point for local label i.
label facePointLabel(const label facei, const label i, const label j) const
Face vertex label offset for a particular i,j,k position.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:50
bool flatFaceOrEdge(const label i, const label j, const label k) const
Return true if point i,j,k addresses a block flat face or edge.
#define w0
Definition: blockCreate.C:30
const Cmpt & x() const
Definition: VectorI.H:75
void correctFacePoints(FixedList< pointField, 6 > &) const
Correct the location of the given face-points.
label nCurvedFaces() const
Number of curved faces in this block.
#define w3
Definition: blockCreate.C:33
label nPoints() const
Return the number of points.
#define w8
Definition: blockCreate.C:40
void setSize(const label)
Reset size of List.
Definition: List.C:281
label patchi
vector point
Point is a vector.
Definition: point.H:41
#define w6
Definition: blockCreate.C:37
#define w10
Definition: blockCreate.C:42
label pointLabel(const label i, const label j, const label k) const
Vertex label offset for a particular i,j,k position.
#define w5
Definition: blockCreate.C:36
#define w9
Definition: blockCreate.C:41
#define w1
Definition: blockCreate.C:31
volScalarField & p
#define w7
Definition: blockCreate.C:38