Alsvinn  0.5.3
The fast FVM simulator with UQ support
volume_foreach.hpp
Go to the documentation of this file.
1 /* Copyright (c) 2018 ETH Zurich, Kjetil Olsen Lye
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation, either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  */
15 
16 #pragma once
17 #include "alsfvm/volume/Volume.hpp"
18 #include <functional>
19 #include <array>
21 #include "alsfvm/grid/Grid.hpp"
25 
26 namespace alsfvm {
27 namespace volume {
28 
29 
43 template<class Function>
44 inline void for_each_cell_index(const Volume& in, const Function& function,
45  ivec3 offsetStart = {0, 0, 0},
46  ivec3 offsetEnd = {0, 0, 0}) {
47  const size_t nx = in.getTotalNumberOfXCells();
48  const size_t ny = in.getTotalNumberOfYCells();
49 
50  const size_t endx = in.getTotalNumberOfXCells() - offsetEnd[0];
51  const size_t endy = in.getTotalNumberOfYCells() - offsetEnd[1];
52  const size_t endz = in.getTotalNumberOfZCells() - offsetEnd[2];
53 
54  for (size_t k = offsetStart[2]; k < endz; k++) {
55  for (size_t j = offsetStart[1]; j < endy; j++) {
56  for (size_t i = offsetStart[0]; i < endx; i++) {
57  size_t index = k * nx * ny + j * nx + i;
58  function(index);
59  }
60  }
61  }
62 }
63 
71 template<size_t direction, bool parallel = false>
73  const std::function<void(size_t leftIndex, size_t middleIndex, size_t rightIndex)>&
74  function,
75  ivec3 offsetStart = {0, 0, 0},
76  ivec3 offsetEnd = {0, 0, 0}) {
77  static_assert(direction < 3, "Direction can be either 0, 1 or 2");
78  const bool xDir = direction == 0;
79  const bool yDir = direction == 1;
80  const bool zDir = direction == 2;
81 
82  const auto view = in.getScalarMemoryArea(0)->getView();
83 
84  const size_t nx = in.getTotalNumberOfXCells() - offsetEnd[0];
85  const size_t ny = in.getTotalNumberOfYCells() - offsetEnd[1];
86  const size_t nz = in.getTotalNumberOfZCells() - offsetEnd[2];
87 
88  for (size_t z = offsetStart[2]; z < nz; z++) {
89  for (size_t y = offsetStart[1]; y < ny; y++) {
90  for (size_t x = offsetStart[0]; x < nx; x++) {
91  const size_t index = view.index(x, y, z);
92  const size_t leftIndex = view.index(x - xDir, y - yDir, z - zDir);
93 
94  const size_t rightIndex = view.index(x + xDir, y + yDir, z + zDir);
95  function(leftIndex, index, rightIndex);
96  }
97  }
98  }
99 }
100 
101 
111 inline void for_each_cell_index_with_neighbours(size_t direction,
112  const Volume& in,
113  const std::function<void(size_t leftIndex, size_t middleIndex, size_t rightIndex)>&
114  function,
115  ivec3 offsetStart = { 0, 0, 0 },
116  ivec3 offsetEnd = { 0, 0, 0 }) {
117  if (direction == 0) {
118  for_each_cell_index_with_neighbours<0>(in, function, offsetStart, offsetEnd);
119  } else if (direction == 1) {
120  for_each_cell_index_with_neighbours<1>(in, function, offsetStart, offsetEnd);
121  } else if (direction == 2) {
122  for_each_cell_index_with_neighbours<2>(in, function, offsetStart, offsetEnd);
123  } else {
124  THROW("Unsupported direction: " << direction);
125  }
126 }
127 
128 
129 
130 
131 template<class VariableStruct>
132 inline VariableStruct expandVariableStruct(const std::array<const real*, 0>& in,
133  size_t index) {
134  // Yes, we want an empty one, since this is called when we do not
135  // have extra variables (eg. burgers)
136  return VariableStruct();
137 }
138 
139 template<class VariableStruct>
140 inline VariableStruct expandVariableStruct(const std::array<const real*, 1>& in,
141  size_t index) {
142  return VariableStruct(in[0][index]);
143 }
144 
145 
146 template<class VariableStruct>
147 inline VariableStruct expandVariableStruct(const std::array<const real*, 3>& in,
148  size_t index) {
149  return VariableStruct(in[0][index], in[1][index], in[2][index]);
150 }
151 
152 
153 template<class VariableStruct>
154 inline VariableStruct expandVariableStruct(const std::array<const real*, 2>& in,
155  size_t index) {
156  return VariableStruct(in[0][index], in[1][index]);
157 }
158 
159 template<class VariableStruct>
160 inline VariableStruct expandVariableStruct(const std::array<const real*, 5>& in,
161  size_t index) {
162  return VariableStruct(in[0][index], in[1][index], in[2][index], in[3][index],
163  in[4][index]);
164 }
165 
166 template<class VariableStruct>
167 inline VariableStruct expandVariableStruct(const std::array<const real*, 4>& in,
168  size_t index) {
169  return VariableStruct(in[0][index], in[1][index], in[2][index], in[3][index]);
170 }
171 
172 
173 template<class VariableStruct>
174 inline void saveVariableStruct(const VariableStruct& in, size_t index,
175  std::array<real*, 1>& out) {
176  real* inAsRealPointer = (real*)&in;
177 
178  out[0][index] = inAsRealPointer[0];
179 
180 }
181 
182 template<class VariableStruct>
183 inline void saveVariableStruct(const VariableStruct& in, size_t index,
184  std::array<real*, 0>& out) {
185  // Yes, we want an empty one, since this is called when we do not
186  // have extra variables (eg. burgers)
187 }
188 
189 
190 template<class VariableStruct>
191 inline void saveVariableStruct(const VariableStruct& in, size_t index,
192  std::array<real*, 4>& out) {
193  real* inAsRealPointer = (real*)&in;
194 
195  out[0][index] = inAsRealPointer[0];
196  out[1][index] = inAsRealPointer[1];
197  out[2][index] = inAsRealPointer[2];
198  out[3][index] = inAsRealPointer[3];
199 }
200 
201 
202 template<class VariableStruct>
203 inline void saveVariableStruct(const VariableStruct& in, size_t index,
204  std::array<real*, 3>& out) {
205  real* inAsRealPointer = (real*)&in;
206 
207  out[0][index] = inAsRealPointer[0];
208  out[1][index] = inAsRealPointer[1];
209  out[2][index] = inAsRealPointer[2];
210 
211 }
212 
213 template<class VariableStruct>
214 inline void saveVariableStruct(const VariableStruct& in, size_t index,
215  std::array<real*, 2>& out) {
216  real* inAsRealPointer = (real*)&in;
217 
218  out[0][index] = inAsRealPointer[0];
219  out[1][index] = inAsRealPointer[1];
220 
221 }
222 
223 
224 
225 template<class VariableStruct>
226 inline void saveVariableStruct(const VariableStruct& in, size_t index,
227  std::array<real*, 5>& out) {
228  real* inAsRealPointer = (real*)&in;
229 
230  out[0][index] = inAsRealPointer[0];
231  out[1][index] = inAsRealPointer[1];
232  out[2][index] = inAsRealPointer[2];
233  out[3][index] = inAsRealPointer[3];
234  out[4][index] = inAsRealPointer[4];
235 }
236 
237 
248 template<class VariableStructIn, class VariableStructOut>
249 inline void transform_volume(const Volume& in, Volume& out,
250  const std::function<VariableStructOut(const VariableStructIn&)>& function) {
251 
252  std::array < const real*, sizeof(VariableStructIn) / sizeof(real) > pointersIn;
253  pointersIn.fill(nullptr);
254 
255  if (pointersIn.size() != in.getNumberOfVariables()) {
256  THROW("Number of variables given does not match struct size."
257  << " Number of variables given: "
258  << in.getNumberOfVariables());
259  }
260 
261  for (size_t i = 0; i < in.getNumberOfVariables(); i++) {
262  pointersIn[i] = in.getScalarMemoryArea(i)->getPointer();
263  }
264 
265  std::array < real*, sizeof(VariableStructOut) / sizeof(real) > pointersOut;
266  pointersOut.fill(nullptr);
267 
268  for (size_t i = 0; i < out.getNumberOfVariables(); i++) {
269  pointersOut[i] = out.getScalarMemoryArea(i)->getPointer();
270  }
271 
272  for_each_cell_index(in, [&](size_t index) {
273  auto out = function(expandVariableStruct<VariableStructIn>(pointersIn, index));
274  saveVariableStruct(out, index, pointersOut);
275 
276  });
277 
278 }
279 
280 
290 template<class VariableStruct>
291 inline void for_each_cell(const Volume& in,
292  const std::function<void(const VariableStruct&, size_t index)>& function) {
293  std::array < const real*, sizeof(VariableStruct) / sizeof(real) > pointersIn;
294  pointersIn.fill(nullptr);
295 
296  if (pointersIn.size() != in.getNumberOfVariables()) {
297  THROW("We expected to get " << pointersIn.size() << " variables, but got " <<
298  in.getNumberOfVariables());
299  }
300 
301  for (size_t i = 0; i < in.getNumberOfVariables(); i++) {
302  pointersIn[i] = in.getScalarMemoryArea(i)->getPointer();
303  }
304 
305  for_each_cell_index(in, [&](size_t index) {
306  function(expandVariableStruct<VariableStruct>(pointersIn, index), index);
307  });
308 }
309 
320 template<class VariableStructA, class VariableStructB>
321 inline void for_each_cell(const Volume& inA, const Volume& inB,
322  const std::function<void(const VariableStructA&, const VariableStructB&, size_t index)>&
323  function) {
324  std::array < const real*, sizeof(VariableStructA) / sizeof(real) > pointersInA;
325  pointersInA.fill(nullptr);
326  std::array < const real*, sizeof(VariableStructB) / sizeof(real) > pointersInB;
327  pointersInB.fill(nullptr);
328 
329  if (pointersInA.size() != inA.getNumberOfVariables()) {
330  THROW("We expected to get " << pointersInA.size() << " variables, but got " <<
331  inA.getNumberOfVariables());
332  }
333 
334  if (pointersInB.size() != inB.getNumberOfVariables()) {
335  THROW("We expected to get " << pointersInB.size() << " variables, but got " <<
336  inB.getNumberOfVariables());
337  }
338 
339 
340  for (size_t i = 0; i < inA.getNumberOfVariables(); i++) {
341  pointersInA[i] = inA.getScalarMemoryArea(i)->getPointer();
342  }
343 
344  for (size_t i = 0; i < inB.getNumberOfVariables(); i++) {
345  pointersInB[i] = inB.getScalarMemoryArea(i)->getPointer();
346  }
347 
348  for_each_cell_index(inA, [&](size_t index) {
349  function(expandVariableStruct<VariableStructA>(pointersInA, index),
350  expandVariableStruct<VariableStructB>(pointersInB, index), index);
351  });
352 }
353 
373 template<class VariableStructA, class VariableStructB>
374 inline void fill_volume(Volume& outA, Volume& outB,
375  const grid::Grid& grid,
376  const std::function < void(real x, real y, real z, VariableStructA& outA,
377  VariableStructB& outB) >& fillerFunction) {
378 
379  std::array < real*, sizeof(VariableStructA) / sizeof(real) > pointersOutA;
380  pointersOutA.fill(nullptr);
381  std::array < real*, sizeof(VariableStructB) / sizeof(real) > pointersOutB;
382  pointersOutB.fill(nullptr);
383 
384  if (pointersOutA.size() != outA.getNumberOfVariables()) {
385  THROW("We expected to get " << pointersOutA.size() << " variables, but got " <<
386  outA.getNumberOfVariables());
387  }
388 
389  if (pointersOutB.size() != outB.getNumberOfVariables()) {
390  THROW("We expected to get " << pointersOutB.size() << " variables, but got " <<
391  outB.getNumberOfVariables());
392  }
393 
394 
395  for (size_t i = 0; i < outA.getNumberOfVariables(); i++) {
396  pointersOutA[i] = outA.getScalarMemoryArea(i)->getPointer();
397  }
398 
399  for (size_t i = 0; i < outB.getNumberOfVariables(); i++) {
400  pointersOutB[i] = outB.getScalarMemoryArea(i)->getPointer();
401  }
402 
403  const size_t nx = outA.getTotalNumberOfXCells();
404  const size_t ny = outA.getTotalNumberOfYCells();
405  const size_t nz = outA.getTotalNumberOfZCells();
406 
407  const size_t nxGrid = grid.getDimensions().x;
408  const size_t nyGrid = grid.getDimensions().y;
409 
410  auto& midPoints = grid.getCellMidpoints();
411 
412  const size_t ghostX = outA.getNumberOfXGhostCells();
413  const size_t ghostY = outA.getNumberOfYGhostCells();
414  const size_t ghostZ = outA.getNumberOfZGhostCells();
415 
416  for (size_t k = ghostZ; k < nz - ghostZ; k++) {
417  for (size_t j = ghostY; j < ny - ghostY; j++) {
418  for (size_t i = ghostX; i < nx - ghostX; i++) {
419  size_t index = k * nx * ny + j * nx + i;
420  size_t midpointIndex = (k - ghostZ) * nxGrid * nyGrid
421  + (j - ghostY) * nxGrid + (i - ghostX);
422  auto midPoint = midPoints[midpointIndex];
423  VariableStructA a;
424  VariableStructB b;
425  fillerFunction(midPoint.x, midPoint.y, midPoint.z, a, b);
426  saveVariableStruct(a, index, pointersOutA);
427  saveVariableStruct(b, index, pointersOutB);
428  }
429  }
430  }
431 
432 
433 }
434 
435 
436 
441 inline void for_each_midpoint(const Volume& volume,
442  const grid::Grid& grid,
443  const std::function
444  < void(real x, real y, real z, size_t index, size_t i, size_t j, size_t k) >&
445  function) {
446 
447  const size_t nx = volume.getTotalNumberOfXCells();
448  const size_t ny = volume.getTotalNumberOfYCells();
449  const size_t nz = volume.getTotalNumberOfZCells();
450 
451  const size_t nxGrid = grid.getDimensions().x;
452  const size_t nyGrid = grid.getDimensions().y;
453 
454  auto& midPoints = grid.getCellMidpoints();
455 
456  const size_t ghostX = volume.getNumberOfXGhostCells();
457  const size_t ghostY = volume.getNumberOfYGhostCells();
458  const size_t ghostZ = volume.getNumberOfZGhostCells();
459 
460  for (size_t k = ghostZ; k < nz - ghostZ; k++) {
461  for (size_t j = ghostY; j < ny - ghostY; j++) {
462  for (size_t i = ghostX; i < nx - ghostX; i++) {
463  size_t index = k * nx * ny + j * nx + i;
464 
465  size_t midpointIndex = (k - ghostZ) * nxGrid * nyGrid
466  + (j - ghostY) * nxGrid + (i - ghostX);
467  auto midPoint = midPoints[midpointIndex];
468 
469  function(midPoint.x, midPoint.y, midPoint.z, index, i - ghostX, j - ghostY,
470  k - ghostZ);
471  }
472  }
473  }
474 
475 
476 }
477 
484 inline void for_each_midpoint(const Volume& volume,
485  const grid::Grid& grid,
486  const std::function < void(real x, real y, real z, size_t index) >& function) {
487  for_each_midpoint(volume, grid, [&](real x, real y, real z, size_t index,
488  size_t, size_t, size_t) {
489  function(x, y, z, index);
490  });
491 }
492 
493 
512 template<class VariableStruct>
513 inline void fill_volume(Volume& out,
514  const grid::Grid& grid,
515  const std::function < void(real x, real y, real z, VariableStruct& out) >&
516  fillerFunction) {
517 
518  std::array < real*, sizeof(VariableStruct) / sizeof(real) > pointersOut;
519  pointersOut.fill(nullptr);
520 
521  if (pointersOut.size() != out.getNumberOfVariables()) {
522  THROW("We expected to get " << pointersOut.size() << " variables, but got " <<
523  out.getNumberOfVariables());
524  }
525 
526 
527 
528  for (size_t i = 0; i < out.getNumberOfVariables(); i++) {
529  pointersOut[i] = out.getScalarMemoryArea(i)->getPointer();
530  }
531 
532  for_each_midpoint(out, grid, [&](real x, real y, real z, size_t index) {
533  VariableStruct a;
534  fillerFunction(x, y, z, a);
535  saveVariableStruct(a, index, pointersOut);
536  });
537 
538 
539 }
540 
541 
542 
564 template<size_t direction>
565 inline void for_each_internal_volume_index(const Volume& volume,
566  const std::function<void(size_t indexLeft, size_t indexMiddle, size_t indexRight)>&
567  function,
568  size_t ghostLayer) {
569  static_assert(direction < 3, "We only support up to three dimensions.");
570  const bool zDir = direction == 2;
571  const bool yDir = direction == 1;
572  const bool xDir = direction == 0;
573 
574  const size_t nx = volume.getTotalNumberOfXCells();
575  const size_t ny = volume.getTotalNumberOfYCells();
576  const size_t nz = volume.getTotalNumberOfZCells();
577 
578  const size_t startZ = ghostLayer * zDir;
579  const size_t startY = ghostLayer * yDir;
580  const size_t startX = ghostLayer * xDir;
581 
582  const size_t endZ = nz - ghostLayer * zDir;
583  const size_t endY = ny - ghostLayer * yDir;
584  const size_t endX = nx - ghostLayer * xDir;
585 
586 
587 
588  for (size_t z = startZ; z < endZ; z++) {
589  for (size_t y = startY; y < endY; y++) {
590  for (size_t x = startX; x < endX; x++) {
591  const size_t index = z * nx * ny + y * nx + x;
592  const size_t leftIndex = (z - zDir) * nx * ny
593  + (y - yDir) * nx
594  + (x - xDir);
595 
596  const size_t rightIndex = (z + zDir) * nx * ny
597  + (y + yDir) * nx
598  + (x + xDir);
599 
600  function(leftIndex, index, rightIndex);
601  }
602  }
603  }
604 }
605 
606 
617 inline void for_each_internal_volume_index(const Volume& volume,
618  const std::function<void( size_t indexMiddle)>& function) {
619 
620 
621  const size_t ngx = volume.getNumberOfXGhostCells();
622  const size_t ngy = volume.getNumberOfYGhostCells();
623  const size_t ngz = volume.getNumberOfZGhostCells();
624  const size_t nx = volume.getTotalNumberOfXCells();
625  const size_t ny = volume.getTotalNumberOfYCells();
626  const size_t nz = volume.getTotalNumberOfZCells();
627 
628  const size_t startZ = ngz;
629  const size_t startY = ngy;
630  const size_t startX = ngx;
631 
632  const size_t endZ = nz - ngz;
633  const size_t endY = ny - ngy;
634  const size_t endX = nx - ngx;
635 
636 
637 
638  for (size_t z = startZ; z < endZ; z++) {
639  for (size_t y = startY; y < endY; y++) {
640  for (size_t x = startX; x < endX; x++) {
641  const size_t index = z * nx * ny + y * nx + x;
642 
643 
644  function(index);
645  }
646  }
647  }
648 }
649 
650 
672 template<size_t direction>
673 inline void for_each_internal_volume_index(const Volume& volume,
674  const std::function<void(size_t indexLeft, size_t indexMiddle, size_t indexRight)>&
675  function) {
676  for_each_internal_volume_index<direction>(volume, function,
677  volume.getNumberOfXGhostCells());
678 }
679 
683 inline void for_each_internal_volume_index(const Volume& volume,
684  size_t direction,
685  const std::function<void(size_t indexLeft, size_t indexMiddle, size_t indexRight)>&
686  function,
687  size_t ghostLayer) {
688 
689  if (direction == 0) {
690  for_each_internal_volume_index<0>(volume, function, ghostLayer);
691  } else if (direction == 1) {
692  for_each_internal_volume_index<1>(volume, function, ghostLayer);
693  } else if (direction == 2) {
694  for_each_internal_volume_index<2>(volume, function, ghostLayer);
695  } else {
696  THROW("We only support direction 0, 1 or 2, was given: " << direction);
697  }
698 }
699 
700 inline void for_each_internal_volume_index(const Volume& volume,
701  size_t direction,
702  const std::function<void(size_t indexLeft, size_t indexMiddle, size_t indexRight)>&
703  function) {
704  for_each_internal_volume_index(volume, direction, function,
705  volume.getNumberOfXGhostCells());
706 }
707 }
708 }
ivec3 getDimensions() const
Definition: Grid.cpp:127
Definition: Grid.hpp:27
#define THROW(message)
Definition: Exception.hpp:27
const std::vector< rvec3 > & getCellMidpoints() const
Definition: Grid.cpp:148
void for_each_cell_index_with_neighbours(const Volume &in, const std::function< void(size_t leftIndex, size_t middleIndex, size_t rightIndex)> &function, ivec3 offsetStart={0, 0, 0}, ivec3 offsetEnd={0, 0, 0})
Definition: volume_foreach.hpp:72
virtual alsfvm::shared_ptr< memory::Memory< real > > & getScalarMemoryArea(size_t index)
getScalarMemoryArea gets the scalar memory area (real)
Definition: Volume.cpp:107
The Volume class represents a volume (a collection of cells with values for each cell (eg...
Definition: Volume.hpp:30
void transform_volume(const Volume &in, Volume &out, const std::function< VariableStructOut(const VariableStructIn &)> &function)
Definition: volume_foreach.hpp:249
virtual size_t getNumberOfVariables() const
getNumberOfVariables gets the number of variables used
Definition: Volume.cpp:83
Definition: vec3.hpp:25
virtual size_t getTotalNumberOfYCells() const
Definition: Volume.cpp:268
VariableStruct expandVariableStruct(const std::array< const real *, 0 > &in, size_t index)
Definition: volume_foreach.hpp:132
size_t nx
Definition: VolumeFactory.cpp:87
double real
Definition: types.hpp:65
virtual size_t getTotalNumberOfZCells() const
Definition: Volume.cpp:275
void saveVariableStruct(const VariableStruct &in, size_t index, std::array< real *, 1 > &out)
Definition: volume_foreach.hpp:174
virtual size_t getNumberOfXGhostCells() const
Definition: Volume.cpp:238
virtual size_t getNumberOfZGhostCells() const
Definition: Volume.cpp:254
void fill_volume(Volume &outA, Volume &outB, const grid::Grid &grid, const std::function< void(real x, real y, real z, VariableStructA &outA, VariableStructB &outB) > &fillerFunction)
Definition: volume_foreach.hpp:374
size_t ny
Definition: VolumeFactory.cpp:88
const grid::Grid & grid
Definition: NumericalFluxFactory.cpp:104
virtual size_t getNumberOfYGhostCells() const
Definition: Volume.cpp:246
#define static_assert(x, y)
Definition: types.hpp:52
void for_each_midpoint(const Volume &volume, const grid::Grid &grid, const std::function< void(real x, real y, real z, size_t index, size_t i, size_t j, size_t k) > &function)
Definition: volume_foreach.hpp:441
size_t nz
Definition: VolumeFactory.cpp:89
T y
Definition: vec3.hpp:27
void for_each_internal_volume_index(const Volume &volume, const std::function< void(size_t indexLeft, size_t indexMiddle, size_t indexRight)> &function, size_t ghostLayer)
Definition: volume_foreach.hpp:565
Various utility functions to implement the tecno flux.
Definition: types.hpp:30
virtual size_t getTotalNumberOfXCells() const
Definition: Volume.cpp:261
T x
Definition: vec3.hpp:26
void for_each_cell(const Volume &in, const std::function< void(const VariableStruct &, size_t index)> &function)
Definition: volume_foreach.hpp:291
void for_each_cell_index(const Volume &in, const Function &function, ivec3 offsetStart={0, 0, 0}, ivec3 offsetEnd={0, 0, 0})
Definition: volume_foreach.hpp:44