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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2016 OpenFOAM Foundation
9  Copyright (C) 2019 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "block.H"
30 
31 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
32 
33 #define w0 w[0][i]
34 #define w1 w[1][i]
35 #define w2 w[2][i]
36 #define w3 w[3][i]
37 
38 #define w4 w[4][j]
39 #define w5 w[5][j]
40 #define w6 w[6][j]
41 #define w7 w[7][j]
42 
43 #define w8 w[8][k]
44 #define w9 w[9][k]
45 #define w10 w[10][k]
46 #define w11 w[11][k]
47 
48 void Foam::block::createPoints()
49 {
50  // Set local variables for mesh specification
51  const label ni = density().x();
52  const label nj = density().y();
53  const label nk = density().z();
54 
55  const point& p000 = blockPoint(0);
56  const point& p100 = blockPoint(1);
57  const point& p110 = blockPoint(2);
58  const point& p010 = blockPoint(3);
59 
60  const point& p001 = blockPoint(4);
61  const point& p101 = blockPoint(5);
62  const point& p111 = blockPoint(6);
63  const point& p011 = blockPoint(7);
64 
65  // List of edge point and weighting factors
66  pointField p[12];
67  scalarList w[12];
68  label nCurvedEdges = edgesPointsWeights(p, w);
69 
70  points_.resize(nPoints());
71 
72  points_[pointLabel(0, 0, 0)] = p000;
73  points_[pointLabel(ni, 0, 0)] = p100;
74  points_[pointLabel(ni, nj, 0)] = p110;
75  points_[pointLabel(0, nj, 0)] = p010;
76  points_[pointLabel(0, 0, nk)] = p001;
77  points_[pointLabel(ni, 0, nk)] = p101;
78  points_[pointLabel(ni, nj, nk)] = p111;
79  points_[pointLabel(0, nj, nk)] = p011;
80 
81  for (label k=0; k<=nk; k++)
82  {
83  for (label j=0; j<=nj; j++)
84  {
85  for (label i=0; i<=ni; i++)
86  {
87  // Skip block vertices
88  if (vertex(i, j, k)) continue;
89 
90  const label vijk = pointLabel(i, j, k);
91 
92  // Calculate the weighting factors for all edges
93 
94  // x-direction
95  scalar wx1 = (1 - w0)*(1 - w4)*(1 - w8) + w0*(1 - w5)*(1 - w9);
96  scalar wx2 = (1 - w1)*w4*(1 - w11) + w1*w5*(1 - w10);
97  scalar wx3 = (1 - w2)*w7*w11 + w2*w6*w10;
98  scalar wx4 = (1 - w3)*(1 - w7)*w8 + w3*(1 - w6)*w9;
99 
100  const scalar sumWx = wx1 + wx2 + wx3 + wx4;
101  wx1 /= sumWx;
102  wx2 /= sumWx;
103  wx3 /= sumWx;
104  wx4 /= sumWx;
105 
106 
107  // y-direction
108  scalar wy1 = (1 - w4)*(1 - w0)*(1 - w8) + w4*(1 - w1)*(1 - w11);
109  scalar wy2 = (1 - w5)*w0*(1 - w9) + w5*w1*(1 - w10);
110  scalar wy3 = (1 - w6)*w3*w9 + w6*w2*w10;
111  scalar wy4 = (1 - w7)*(1 - w3)*w8 + w7*(1 - w2)*w11;
112 
113  const scalar sumWy = wy1 + wy2 + wy3 + wy4;
114  wy1 /= sumWy;
115  wy2 /= sumWy;
116  wy3 /= sumWy;
117  wy4 /= sumWy;
118 
119 
120  // z-direction
121  scalar wz1 = (1 - w8)*(1 - w0)*(1 - w4) + w8*(1 - w3)*(1 - w7);
122  scalar wz2 = (1 - w9)*w0*(1 - w5) + w9*w3*(1 - w6);
123  scalar wz3 = (1 - w10)*w1*w5 + w10*w2*w6;
124  scalar wz4 = (1 - w11)*(1 - w1)*w4 + w11*(1 - w2)*w7;
125 
126  const scalar sumWz = wz1 + wz2 + wz3 + wz4;
127  wz1 /= sumWz;
128  wz2 /= sumWz;
129  wz3 /= sumWz;
130  wz4 /= sumWz;
131 
132 
133  // Points on straight edges
134  const vector edgex1 = p000 + (p100 - p000)*w0;
135  const vector edgex2 = p010 + (p110 - p010)*w1;
136  const vector edgex3 = p011 + (p111 - p011)*w2;
137  const vector edgex4 = p001 + (p101 - p001)*w3;
138 
139  const vector edgey1 = p000 + (p010 - p000)*w4;
140  const vector edgey2 = p100 + (p110 - p100)*w5;
141  const vector edgey3 = p101 + (p111 - p101)*w6;
142  const vector edgey4 = p001 + (p011 - p001)*w7;
143 
144  const vector edgez1 = p000 + (p001 - p000)*w8;
145  const vector edgez2 = p100 + (p101 - p100)*w9;
146  const vector edgez3 = p110 + (p111 - p110)*w10;
147  const vector edgez4 = p010 + (p011 - p010)*w11;
148 
149  // Add the contributions
150  points_[vijk] =
151  (
152  wx1*edgex1 + wx2*edgex2 + wx3*edgex3 + wx4*edgex4
153  + wy1*edgey1 + wy2*edgey2 + wy3*edgey3 + wy4*edgey4
154  + wz1*edgez1 + wz2*edgez2 + wz3*edgez3 + wz4*edgez4
155  )/3;
156 
157 
158  // Apply curved-edge correction if block has curved edges
159  if (nCurvedEdges)
160  {
161  // Calculate the correction vectors
162  const vector corx1 = wx1*(p[0][i] - edgex1);
163  const vector corx2 = wx2*(p[1][i] - edgex2);
164  const vector corx3 = wx3*(p[2][i] - edgex3);
165  const vector corx4 = wx4*(p[3][i] - edgex4);
166 
167  const vector cory1 = wy1*(p[4][j] - edgey1);
168  const vector cory2 = wy2*(p[5][j] - edgey2);
169  const vector cory3 = wy3*(p[6][j] - edgey3);
170  const vector cory4 = wy4*(p[7][j] - edgey4);
171 
172  const vector corz1 = wz1*(p[8][k] - edgez1);
173  const vector corz2 = wz2*(p[9][k] - edgez2);
174  const vector corz3 = wz3*(p[10][k] - edgez3);
175  const vector corz4 = wz4*(p[11][k] - edgez4);
176 
177  points_[vijk] +=
178  (
179  corx1 + corx2 + corx3 + corx4
180  + cory1 + cory2 + cory3 + cory4
181  + corz1 + corz2 + corz3 + corz4
182  );
183  }
184  }
185  }
186  }
187 
188  if (!nCurvedFaces()) return;
189 
190  // Apply curvature correction to face points
191  FixedList<pointField, 6> facePoints(this->facePoints(points_));
193 
194  // Loop over points and apply the correction from the i-faces
195  for (label ii=0; ii<=ni; ii++)
196  {
197  // Process the points on the curved faces last
198  const label i = (ii + 1)%(ni + 1);
199 
200  for (label j=0; j<=nj; j++)
201  {
202  for (label k=0; k<=nk; k++)
203  {
204  // Skip non-curved faces and edges
205  if (flatFaceOrEdge(i, j, k)) continue;
206 
207  const label vijk = pointLabel(i, j, k);
208 
209  scalar wf0 =
210  (
211  (1 - w0)*(1 - w4)*(1 - w8)
212  + (1 - w1)*w4*(1 - w11)
213  + (1 - w2)*w7*w11
214  + (1 - w3)*(1 - w7)*w8
215  );
216 
217  scalar wf1 =
218  (
219  w0*(1 - w5)*(1 - w9)
220  + w1*w5*(1 - w10)
221  + w2*w5*w10
222  + w3*(1 - w6)*w9
223  );
224 
225  const scalar sumWf = wf0 + wf1;
226  wf0 /= sumWf;
227  wf1 /= sumWf;
228 
229  points_[vijk] +=
230  (
231  wf0
232  *(
233  facePoints[0][facePointLabel(0, j, k)]
234  - points_[pointLabel(0, j, k)]
235  )
236  + wf1
237  *(
238  facePoints[1][facePointLabel(1, j, k)]
239  - points_[pointLabel(ni, j, k)]
240  )
241  );
242  }
243  }
244  }
245 
246  // Loop over points and apply the correction from the j-faces
247  for (label jj=0; jj<=nj; jj++)
248  {
249  // Process the points on the curved faces last
250  const label j = (jj + 1)%(nj + 1);
251 
252  for (label i=0; i<=ni; i++)
253  {
254  for (label k=0; k<=nk; k++)
255  {
256  // Skip flat faces and edges
257  if (flatFaceOrEdge(i, j, k)) continue;
258 
259  const label vijk = pointLabel(i, j, k);
260 
261  scalar wf2 =
262  (
263  (1 - w4)*(1 - w1)*(1 - w8)
264  + (1 - w5)*w0*(1 - w9)
265  + (1 - w6)*w3*w9
266  + (1 - w7)*(1 - w3)*w8
267  );
268 
269  scalar wf3 =
270  (
271  w4*(1 - w1)*(1 - w11)
272  + w5*w1*(1 - w10)
273  + w6*w2*w10
274  + w7*(1 - w2)*w11
275  );
276 
277  const scalar sumWf = wf2 + wf3;
278  wf2 /= sumWf;
279  wf3 /= sumWf;
280 
281  points_[vijk] +=
282  (
283  wf2
284  *(
285  facePoints[2][facePointLabel(2, i, k)]
286  - points_[pointLabel(i, 0, k)]
287  )
288  + wf3
289  *(
290  facePoints[3][facePointLabel(3, i, k)]
291  - points_[pointLabel(i, nj, k)]
292  )
293  );
294  }
295  }
296  }
297 
298  // Loop over points and apply the correction from the k-faces
299  for (label kk=0; kk<=nk; kk++)
300  {
301  // Process the points on the curved faces last
302  const label k = (kk + 1)%(nk + 1);
303 
304  for (label i=0; i<=ni; i++)
305  {
306  for (label j=0; j<=nj; j++)
307  {
308  // Skip flat faces and edges
309  if (flatFaceOrEdge(i, j, k)) continue;
310 
311  const label vijk = pointLabel(i, j, k);
312 
313  scalar wf4 =
314  (
315  (1 - w8)*(1 - w0)*(1 - w4)
316  + (1 - w9)*w0*(1 - w5)
317  + (1 - w10)*w1*w5
318  + (1 - w11)*(1 - w1)*w4
319  );
320 
321  scalar wf5 =
322  (
323  w8*(1 - w3)*(1 - w7)
324  + w9*w3*(1 - w6)
325  + w10*w2*w6
326  + w11*(1 - w2)*w7
327  );
328 
329  const scalar sumWf = wf4 + wf5;
330  wf4 /= sumWf;
331  wf5 /= sumWf;
332 
333  points_[vijk] +=
334  (
335  wf4
336  *(
337  facePoints[4][facePointLabel(4, i, j)]
338  - points_[pointLabel(i, j, 0)]
339  )
340  + wf5
341  *(
342  facePoints[5][facePointLabel(5, i, j)]
343  - points_[pointLabel(i, j, nk)]
344  )
345  );
346  }
347  }
348  }
349 }
350 
351 
352 void Foam::block::createCells()
353 {
354  const label ni = density().x();
355  const label nj = density().y();
356  const label nk = density().z();
357 
358  blockCells_.resize(nCells()); // (ni*nj*nk)
359 
360  label celli = 0;
361 
362  for (label k=0; k<nk; ++k)
363  {
364  for (label j=0; j<nj; ++j)
365  {
366  for (label i=0; i<ni; ++i)
367  {
368  blockCells_[celli][0] = pointLabel(i, j, k);
369  blockCells_[celli][1] = pointLabel(i+1, j, k);
370  blockCells_[celli][2] = pointLabel(i+1, j+1, k);
371  blockCells_[celli][3] = pointLabel(i, j+1, k);
372  blockCells_[celli][4] = pointLabel(i, j, k+1);
373  blockCells_[celli][5] = pointLabel(i+1, j, k+1);
374  blockCells_[celli][6] = pointLabel(i+1, j+1, k+1);
375  blockCells_[celli][7] = pointLabel(i, j+1, k+1);
376 
377  ++celli;
378  }
379  }
380  }
381 }
382 
383 
384 template<class OutputIterator>
385 OutputIterator Foam::block::addBoundaryFaces
386 (
387  const direction shapeFacei,
388  OutputIterator iter
389 ) const
390 {
391  const label ni = density().x();
392  const label nj = density().y();
393  const label nk = density().z();
394 
395  switch (shapeFacei)
396  {
397  // Face 0 == x-min
398  case 0:
399  {
400  for (label k=0; k<nk; ++k)
401  {
402  for (label j=0; j<nj; ++j)
403  {
404  auto& f = *iter;
405  ++iter;
406  f.resize(4);
407 
408  f[0] = pointLabel(0, j, k);
409  f[1] = pointLabel(0, j, k+1);
410  f[2] = pointLabel(0, j+1, k+1);
411  f[3] = pointLabel(0, j+1, k);
412  }
413  }
414  }
415  break;
416 
417  // Face 1 == x-max
418  case 1:
419  {
420  for (label k=0; k<nk; ++k)
421  {
422  for (label j=0; j<nj; ++j)
423  {
424  auto& f = *iter;
425  ++iter;
426  f.resize(4);
427 
428  f[0] = pointLabel(ni, j, k);
429  f[1] = pointLabel(ni, j+1, k);
430  f[2] = pointLabel(ni, j+1, k+1);
431  f[3] = pointLabel(ni, j, k+1);
432  }
433  }
434  }
435  break;
436 
437  // Face 2 == y-min
438  case 2:
439  {
440  for (label i=0; i<ni; ++i)
441  {
442  for (label k=0; k<nk; ++k)
443  {
444  auto& f = *iter;
445  ++iter;
446  f.resize(4);
447 
448  f[0] = pointLabel(i, 0, k);
449  f[1] = pointLabel(i+1, 0, k);
450  f[2] = pointLabel(i+1, 0, k+1);
451  f[3] = pointLabel(i, 0, k+1);
452  }
453  }
454  }
455  break;
456 
457  // Face 3 == y-max
458  case 3:
459  {
460  for (label i=0; i<ni; ++i)
461  {
462  for (label k=0; k<nk; ++k)
463  {
464  auto& f = *iter;
465  ++iter;
466  f.resize(4);
467 
468  f[0] = pointLabel(i, nj, k);
469  f[1] = pointLabel(i, nj, k+1);
470  f[2] = pointLabel(i+1, nj, k+1);
471  f[3] = pointLabel(i+1, nj, k);
472  }
473  }
474  }
475  break;
476 
477  // Face 4 == z-min
478  case 4:
479  {
480  for (label i=0; i<ni; ++i)
481  {
482  for (label j=0; j<nj; ++j)
483  {
484  auto& f = *iter;
485  ++iter;
486  f.resize(4);
487 
488  f[0] = pointLabel(i, j, 0);
489  f[1] = pointLabel(i, j+1, 0);
490  f[2] = pointLabel(i+1, j+1, 0);
491  f[3] = pointLabel(i+1, j, 0);
492  }
493  }
494  }
495  break;
496 
497  // Face 5 == z-max
498  case 5:
499  {
500  for (label i=0; i<ni; ++i)
501  {
502  for (label j=0; j<nj; ++j)
503  {
504  auto& f = *iter;
505  ++iter;
506  f.resize(4);
507 
508  f[0] = pointLabel(i, j, nk);
509  f[1] = pointLabel(i+1, j, nk);
510  f[2] = pointLabel(i+1, j+1, nk);
511  f[3] = pointLabel(i, j+1, nk);
512  }
513  }
514  }
515  break;
516  }
517 
518  return iter;
519 }
520 
521 
522 void Foam::block::createBoundary()
523 {
524  const label countx = (density().y() * density().z());
525  const label county = (density().z() * density().x());
526  const label countz = (density().x() * density().y());
527 
528  direction patchi = 0;
529 
530  // 0 == x-min
531  blockPatches_[patchi].resize(countx);
532  addBoundaryFaces(patchi, blockPatches_[patchi].begin());
533  ++patchi;
534 
535  // 1 == x-max
536  blockPatches_[patchi].resize(countx);
537  addBoundaryFaces(patchi, blockPatches_[patchi].begin());
538  ++patchi;
539 
540  // 2 == y-min
541  blockPatches_[patchi].resize(county);
542  addBoundaryFaces(patchi, blockPatches_[patchi].begin());
543  ++patchi;
544 
545  // 3 == y-max
546  blockPatches_[patchi].resize(county);
547  addBoundaryFaces(patchi, blockPatches_[patchi].begin());
548  ++patchi;
549 
550  // 4 == z-min
551  blockPatches_[patchi].resize(countz);
552  addBoundaryFaces(patchi, blockPatches_[patchi].begin());
553  ++patchi;
554 
555  // 5 == z-max
556  blockPatches_[patchi].resize(countz);
557  addBoundaryFaces(patchi, blockPatches_[patchi].begin());
558  ++patchi;
559 }
560 
561 
562 // ************************************************************************* //
w4
#define w4
Definition: blockCreate.C:38
w11
#define w11
Definition: blockCreate.C:46
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::scalarList
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:64
Foam::Vector::x
const Cmpt & x() const
Access to the vector x component.
Definition: VectorI.H:73
Foam::blockDescriptor::blockPoint
const point & blockPoint(const label i) const
Return block point for local label i.
Definition: blockDescriptorI.H:81
p
volScalarField & p
Definition: createFieldRefs.H:8
stdFoam::begin
constexpr auto begin(C &c) -> decltype(c.begin())
Return iterator to the beginning of the container c.
Definition: stdFoam.H:97
w1
#define w1
Definition: blockCreate.C:34
w7
#define w7
Definition: blockCreate.C:41
Foam::blockDescriptor::facePointLabel
label facePointLabel(const direction facei, const label i, const label j) const
Definition: blockDescriptorI.H:88
Foam::blockDescriptor::vertex
bool vertex(const label i, const label j, const label k) const
Return true if point i,j,k addresses a block vertex.
Definition: blockDescriptorI.H:125
Foam::blockDescriptor::nCurvedFaces
label nCurvedFaces() const
Number of curved faces in this block.
Definition: blockDescriptorI.H:75
w3
#define w3
Definition: blockCreate.C:36
Foam::Vector::z
const Cmpt & z() const
Access to the vector z component.
Definition: VectorI.H:85
w6
#define w6
Definition: blockCreate.C:40
Foam::ijkMesh::pointLabel
label pointLabel(const label i, const label j, const label k) const
The linear point index for an i-j-k position.
Definition: ijkMeshI.H:183
w10
#define w10
Definition: blockCreate.C:45
Foam::blockDescriptor::facePoints
FixedList< pointField, 6 > facePoints(const pointField &points) const
Return the list of face-points for all of the faces of the block.
Definition: blockDescriptor.C:287
Foam::List::resize
void resize(const label newSize)
Adjust allocated size of list.
Definition: ListI.H:139
Foam::blockDescriptor::correctFacePoints
void correctFacePoints(FixedList< pointField, 6 > &) const
Correct the location of the given face-points.
Definition: blockDescriptor.C:343
Foam::blockDescriptor::flatFaceOrEdge
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.
Definition: blockDescriptorI.H:151
w5
#define w5
Definition: blockCreate.C:39
w2
#define w2
Definition: blockCreate.C:35
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::blockDescriptor::edgesPointsWeights
label edgesPointsWeights(pointField(&edgePoints)[12], scalarList(&edgeWeights)[12]) const
Calculate the points and weights for all edges.
Definition: blockDescriptorEdges.C:115
Foam::Vector::y
const Cmpt & y() const
Access to the vector y component.
Definition: VectorI.H:79
f
labelList f(nPoints)
Foam::ijkMesh::nPoints
label nPoints() const
The number of mesh points (nx+1)*(ny+1)*(nz+1) in the i-j-k mesh.
Definition: ijkMeshI.H:55
w8
#define w8
Definition: blockCreate.C:43
k
label k
Boltzmann constant.
Definition: LISASMDCalcMethod2.H:41
Foam::direction
uint8_t direction
Definition: direction.H:47
Foam::blockDescriptor::density
const labelVector & density() const
Return the mesh density (number of cells) in the i,j,k directions.
Definition: blockDescriptorI.H:49
w0
#define w0
Definition: blockCreate.C:33
Foam::point
vector point
Point is a vector.
Definition: point.H:43
block.H
w9
#define w9
Definition: blockCreate.C:44