ThunderEgg  1.0.0
DomainTools.h
Go to the documentation of this file.
1 /***************************************************************************
2  * ThunderEgg, a library for solvers on adaptively refined block-structured
3  * Cartesian grids.
4  *
5  * Copyright (c) 2019-2021 Scott Aiton
6  *
7  * This program is free software: you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation, either version 3 of the License, or
10  * (at your option) any later version.
11  *
12  * This program is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with this program. If not, see <https://www.gnu.org/licenses/>.
19  ***************************************************************************/
20 
21 #ifndef THUNDEREGG_DOMAIN_TOOLS_H
22 #define THUNDEREGG_DOMAIN_TOOLS_H
23 
29 #include <ThunderEgg/RuntimeError.h>
30 #include <ThunderEgg/Vector.h>
31 
32 namespace ThunderEgg {
36 class DomainTools
37 {
38 private:
50  template<int D, typename T>
51  static void _SetValues(const Domain<D>& domain, Vector<D>& vec, int component_index, T func)
52  {
53  SetValues(domain, vec, component_index, func);
54  }
68  template<int D, typename T, typename... Args>
69  static void _SetValues(const Domain<D>& domain,
70  Vector<D>& vec,
71  int component_index,
72  T func,
73  Args... args)
74  {
75  SetValues(domain, vec, component_index, func);
76  _SetValues(domain, vec, component_index + 1, args...);
77  }
89  template<int D, typename T>
90  static void _SetValuesWithGhost(const Domain<D>& domain,
91  Vector<D>& vec,
92  int component_index,
93  T func)
94  {
95  SetValuesWithGhost(domain, vec, component_index, func);
96  }
110  template<int D, typename T, typename... Args>
111  static void _SetValuesWithGhost(const Domain<D>& domain,
112  Vector<D>& vec,
113  int component_index,
114  T func,
115  Args... args)
116  {
117  SetValuesWithGhost(domain, vec, component_index, func);
118  _SetValuesWithGhost(domain, vec, component_index + 1, args...);
119  }
120 
121 public:
133  template<int D>
134  static void GetRealCoord(const PatchInfo<D>& pinfo,
135  const std::array<int, D>& coord,
136  std::array<double, D>& real_coord)
137  {
138  Loop::Unroll<0, D - 1>([&](int dir) {
139  if (coord[dir] == -1) {
140  real_coord[dir] = pinfo.starts[dir];
141  } else if (coord[dir] == pinfo.ns[dir]) {
142  real_coord[dir] = pinfo.starts[dir] + pinfo.spacings[dir] * pinfo.ns[dir];
143  } else {
144  real_coord[dir] =
145  pinfo.starts[dir] + pinfo.spacings[dir] / 2.0 + pinfo.spacings[dir] * coord[dir];
146  }
147  });
148  }
157  template<int D>
158  static void GetRealCoordGhost(const PatchInfo<D>& pinfo,
159  const std::array<int, D>& coord,
160  std::array<double, D>& real_coord)
161  {
162  Loop::Unroll<0, D - 1>([&](int dir) {
163  real_coord[dir] =
164  pinfo.starts[dir] + pinfo.spacings[dir] / 2.0 + pinfo.spacings[dir] * coord[dir];
165  });
166  }
177  template<int D>
178  static void GetRealCoordBound(const PatchInfo<D>& pinfo,
179  const std::array<int, D - 1>& coord,
180  Side<D> s,
181  std::array<double, D>& real_coord)
182  {
183  for (size_t dir = 0; dir < s.getAxisIndex(); dir++) {
184  if (coord[dir] == -1) {
185  real_coord[dir] = pinfo.starts[dir];
186  } else if (coord[dir] == pinfo.ns[dir]) {
187  real_coord[dir] = pinfo.starts[dir] + pinfo.spacings[dir] * pinfo.ns[dir];
188  } else {
189  real_coord[dir] =
190  pinfo.starts[dir] + pinfo.spacings[dir] / 2.0 + pinfo.spacings[dir] * coord[dir];
191  }
192  }
193  if (s.isLowerOnAxis()) {
194  real_coord[s.getAxisIndex()] = pinfo.starts[s.getAxisIndex()];
195  } else {
196  real_coord[s.getAxisIndex()] = pinfo.starts[s.getAxisIndex()] +
197  pinfo.spacings[s.getAxisIndex()] * pinfo.ns[s.getAxisIndex()];
198  }
199  for (size_t dir = s.getAxisIndex() + 1; dir < D; dir++) {
200  if (coord[dir - 1] == -1) {
201  real_coord[dir] = pinfo.starts[dir];
202  } else if (coord[dir - 1] == pinfo.ns[dir]) {
203  real_coord[dir] = pinfo.starts[dir] + pinfo.spacings[dir] * pinfo.ns[dir];
204  } else {
205  real_coord[dir] =
206  pinfo.starts[dir] + pinfo.spacings[dir] / 2.0 + pinfo.spacings[dir] * coord[dir - 1];
207  }
208  }
209  }
219  template<int D>
220  static void SetValues(const Domain<D>& domain,
221  Vector<D>& vec,
222  int component_index,
223  std::function<double(const std::array<double, (int)D>&)> func)
224  {
225  if (component_index >= vec.getNumComponents()) {
226  throw RuntimeError("Invalid component to set");
227  }
228  std::array<double, D> real_coord;
229  for (int i = 0; i < vec.getNumLocalPatches(); i++) {
230  ComponentView<double, D> ld = vec.getComponentView(component_index, i);
231  auto pinfo = domain.getPatchInfoVector()[i];
232  Loop::Nested<D>(ld.getStart(), ld.getEnd(), [&](const std::array<int, D>& coord) {
233  GetRealCoord<D>(pinfo, coord, real_coord);
234  ld[coord] = func(real_coord);
235  });
236  }
237  }
246  static void SetValues(const Domain<3>& domain,
247  Vector<3>& vec,
248  int component_index,
249  std::function<double(double, double, double)> func)
250  {
251  if (component_index >= vec.getNumComponents()) {
252  throw RuntimeError("Invalid component to set");
253  }
254  for (int i = 0; i < vec.getNumLocalPatches(); i++) {
255  ComponentView<double, 3> ld = vec.getComponentView(component_index, i);
256  const PatchInfo<3>& pinfo = domain.getPatchInfoVector()[i];
257  double dx = pinfo.spacings[0];
258  double dy = pinfo.spacings[1];
259  double dz = pinfo.spacings[2];
260  for (int zi = 0; zi < pinfo.ns[2]; zi++) {
261  double z = pinfo.starts[2] + 0.5 * dz + zi * dz;
262  for (int yi = 0; yi < pinfo.ns[1]; yi++) {
263  double y = pinfo.starts[1] + 0.5 * dy + yi * dy;
264  for (int xi = 0; xi < pinfo.ns[0]; xi++) {
265  double x = pinfo.starts[0] + 0.5 * dx + xi * dx;
266  ld(xi, yi, zi) = func(x, y, z);
267  }
268  }
269  }
270  }
271  }
280  static void SetValues(const Domain<2>& domain,
281  Vector<2>& vec,
282  int component_index,
283  std::function<double(double, double)> func)
284  {
285  if (component_index >= vec.getNumComponents()) {
286  throw RuntimeError("Invalid component to set");
287  }
288  for (int i = 0; i < vec.getNumLocalPatches(); i++) {
289  ComponentView<double, 2> ld = vec.getComponentView(component_index, i);
290  const PatchInfo<2>& pinfo = domain.getPatchInfoVector()[i];
291  double dx = pinfo.spacings[0];
292  double dy = pinfo.spacings[1];
293  for (int yi = 0; yi < pinfo.ns[1]; yi++) {
294  double y = pinfo.starts[1] + 0.5 * dy + yi * dy;
295  for (int xi = 0; xi < pinfo.ns[0]; xi++) {
296  double x = pinfo.starts[0] + 0.5 * dx + xi * dx;
297  ld(xi, yi) = func(x, y);
298  }
299  }
300  }
301  }
310  static void SetValues(const Domain<1>& domain,
311  Vector<1>& vec,
312  int component_index,
313  std::function<double(double)> func)
314  {
315  if (component_index >= vec.getNumComponents()) {
316  throw RuntimeError("Invalid component to set");
317  }
318  for (int i = 0; i < vec.getNumLocalPatches(); i++) {
319  ComponentView<double, 1> ld = vec.getComponentView(component_index, i);
320  const PatchInfo<1>& pinfo = domain.getPatchInfoVector()[i];
321  double dx = pinfo.spacings[0];
322  for (int xi = 0; xi < pinfo.ns[0]; xi++) {
323  double x = pinfo.starts[0] + 0.5 * dx + xi * dx;
324  ld(xi) = func(x);
325  }
326  }
327  }
339  template<int D, typename... Args>
340  static void SetValues(const Domain<D>& domain,
341  Vector<D>& vec,
342  std::function<double(const std::array<double, D>&)> func,
343  Args... args)
344  {
345  _SetValues(domain, vec, 0, func, args...);
346  }
356  template<int D>
357  static void SetValuesWithGhost(const Domain<D>& domain,
358  Vector<D>& vec,
359  int component_index,
360  std::function<double(const std::array<double, (int)D>&)> func)
361  {
362  if (component_index >= vec.getNumComponents()) {
363  throw RuntimeError("Invalid component to set");
364  }
365  std::array<double, D> real_coord;
366  for (int i = 0; i < vec.getNumLocalPatches(); i++) {
367  ComponentView<double, D> ld = vec.getComponentView(component_index, i);
368  auto pinfo = domain.getPatchInfoVector()[i];
369  Loop::Nested<D>(ld.getGhostStart(), ld.getGhostEnd(), [&](const std::array<int, D>& coord) {
370  GetRealCoordGhost<D>(pinfo, coord, real_coord);
371  ld[coord] = func(real_coord);
372  });
373  }
374  }
383  static void SetValuesWithGhost(const Domain<3>& domain,
384  Vector<3>& vec,
385  int component_index,
386  std::function<double(double, double, double)> func)
387  {
388  if (component_index >= vec.getNumComponents()) {
389  throw RuntimeError("Invalid component to set");
390  }
391  for (int i = 0; i < vec.getNumLocalPatches(); i++) {
392  ComponentView<double, 3> ld = vec.getComponentView(component_index, i);
393  const PatchInfo<3>& pinfo = domain.getPatchInfoVector()[i];
394  int num_ghost = pinfo.num_ghost_cells;
395  double dx = pinfo.spacings[0];
396  double dy = pinfo.spacings[1];
397  double dz = pinfo.spacings[2];
398  for (int zi = -num_ghost; zi < pinfo.ns[2] + num_ghost; zi++) {
399  double z = pinfo.starts[2] + 0.5 * dz + zi * dz;
400  for (int yi = -num_ghost; yi < pinfo.ns[1] + num_ghost; yi++) {
401  double y = pinfo.starts[1] + 0.5 * dy + yi * dy;
402  for (int xi = -num_ghost; xi < pinfo.ns[0] + num_ghost; xi++) {
403  double x = pinfo.starts[0] + 0.5 * dx + xi * dx;
404  ld(xi, yi, zi) = func(x, y, z);
405  }
406  }
407  }
408  }
409  }
418  static void SetValuesWithGhost(const Domain<2>& domain,
419  Vector<2>& vec,
420  int component_index,
421  std::function<double(double, double)> func)
422  {
423  if (component_index >= vec.getNumComponents()) {
424  throw RuntimeError("Invalid component to set");
425  }
426  for (int i = 0; i < vec.getNumLocalPatches(); i++) {
427  ComponentView<double, 2> ld = vec.getComponentView(component_index, i);
428  const PatchInfo<2>& pinfo = domain.getPatchInfoVector()[i];
429  int num_ghost = pinfo.num_ghost_cells;
430  double dx = pinfo.spacings[0];
431  double dy = pinfo.spacings[1];
432  for (int yi = -num_ghost; yi < pinfo.ns[1] + num_ghost; yi++) {
433  double y = pinfo.starts[1] + 0.5 * dy + yi * dy;
434  for (int xi = -num_ghost; xi < pinfo.ns[0] + num_ghost; xi++) {
435  double x = pinfo.starts[0] + 0.5 * dx + xi * dx;
436  ld(xi, yi) = func(x, y);
437  }
438  }
439  }
440  }
449  static void SetValuesWithGhost(const Domain<1>& domain,
450  Vector<1>& vec,
451  int component_index,
452  std::function<double(double)> func)
453  {
454  if (component_index >= vec.getNumComponents()) {
455  throw RuntimeError("Invalid component to set");
456  }
457  for (int i = 0; i < vec.getNumLocalPatches(); i++) {
458  ComponentView<double, 1> ld = vec.getComponentView(component_index, i);
459  const PatchInfo<1>& pinfo = domain.getPatchInfoVector()[i];
460  int num_ghost = pinfo.num_ghost_cells;
461  double dx = pinfo.spacings[0];
462  for (int xi = -num_ghost; xi < pinfo.ns[0] + num_ghost; xi++) {
463  double x = pinfo.starts[0] + 0.5 * dx + xi * dx;
464  ld(xi) = func(x);
465  }
466  }
467  }
480  template<int D, typename... Args>
481  static void SetValuesWithGhost(const Domain<D>& domain,
482  Vector<D>& vec,
483  std::function<double(const std::array<double, D>&)> func,
484  Args... args)
485  {
486  _SetValuesWithGhost(domain, vec, 0, func, args...);
487  }
498  template<typename... Args>
499  static void SetValuesWithGhost(const Domain<3>& domain,
500  Vector<3>& vec,
501  std::function<double(double, double, double)> func,
502  Args... args)
503  {
504  _SetValuesWithGhost(domain, vec, 0, func, args...);
505  }
512  template<int D>
513  static double Integrate(const Domain<D>& domain, const Vector<D>& u)
514  {
515  double sum = 0;
516 
517  for (const auto& pinfo : domain.getPatchInfoVector()) {
518  for (int c = 0; c < u.getNumComponents(); c++) {
519  ComponentView<const double, D> u_data = u.getComponentView(c, pinfo.local_index);
520 
521  double patch_sum = 0;
522  Loop::Nested<D>(u_data.getStart(), u_data.getEnd(), [&](std::array<int, D> coord) {
523  patch_sum += u_data[coord];
524  });
525 
526  for (size_t i = 0; i < D; i++) {
527  patch_sum *= pinfo.spacings[i];
528  }
529  sum += patch_sum;
530  }
531  }
532  double retval;
533  MPI_Allreduce(&sum, &retval, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
534  return retval;
535  }
536 };
537 } // namespace ThunderEgg
538 #endif
ThunderEgg::Vector::getComponentView
ComponentView< double, D > getComponentView(int component_index, int patch_local_index)
Get the ComponentView for the specified patch and component.
Definition: Vector.h:336
RuntimeError.h
RuntimeError struct.
ThunderEgg::Face::getAxisIndex
auto getAxisIndex() const -> typename std::enable_if< D<=3 &&D >=1 &&M==D - 1 &&N==N, size_t >::type
Get the axis index of this side.
Definition: Face.h:458
ThunderEgg::Loop::Unroll
static void Unroll(T lambda)
Unroll a fixed-length loop.
Definition: Loops.h:72
Vector.h
Vector class.
ThunderEgg::DomainTools::SetValuesWithGhost
static void SetValuesWithGhost(const Domain< 3 > &domain, Vector< 3 > &vec, int component_index, std::function< double(double, double, double)> func)
Set the values (including ghost values) for a vector with the given function.
Definition: DomainTools.h:383
ThunderEgg::View::getGhostStart
const std::array< int, D > & getGhostStart() const
Get the coordinate of the first ghost cell element.
Definition: View.h:215
ThunderEgg::DomainTools::SetValuesWithGhost
static void SetValuesWithGhost(const Domain< D > &domain, Vector< D > &vec, int component_index, std::function< double(const std::array< double,(int) D > &)> func)
Set the values (including ghost values) for a vector with the given function.
Definition: DomainTools.h:357
ThunderEgg::PatchInfo::ns
std::array< int, D > ns
The number of cells in each direction.
Definition: PatchInfo.h:117
ThunderEgg::ComponentView< double, D >
ThunderEgg::Domain
Uses a collection of PatchInfo objects to represent the domain of the problem.
Definition: Domain.h:50
ThunderEgg::DomainTools::SetValuesWithGhost
static void SetValuesWithGhost(const Domain< 2 > &domain, Vector< 2 > &vec, int component_index, std::function< double(double, double)> func)
Set the values (including ghost values) for a vector with the given function.
Definition: DomainTools.h:418
ThunderEgg::Face::isLowerOnAxis
auto isLowerOnAxis() const -> typename std::enable_if< D<=3 &&D >=1 &&M==D - 1 &&N==N, bool >::type
Return if this side is lower on it's axis.
Definition: Face.h:468
ThunderEgg::DomainTools::GetRealCoordGhost
static void GetRealCoordGhost(const PatchInfo< D > &pinfo, const std::array< int, D > &coord, std::array< double, D > &real_coord)
Given a path info object, get the coordinate from a given index into the patch.
Definition: DomainTools.h:158
ThunderEgg
The ThunderEgg namespace.
Definition: BiLinearGhostFiller.h:31
ThunderEgg::DomainTools::SetValues
static void SetValues(const Domain< D > &domain, Vector< D > &vec, std::function< double(const std::array< double, D > &)> func, Args... args)
Set the values for a vector with the given functions.
Definition: DomainTools.h:340
ThunderEgg::View::getEnd
const std::array< int, D > & getEnd() const
Get the coordinate of the last element.
Definition: View.h:211
ThunderEgg::PatchInfo
Contains metadata for a patch.
Definition: PatchInfo.h:51
ThunderEgg::DomainTools::GetRealCoord
static void GetRealCoord(const PatchInfo< D > &pinfo, const std::array< int, D > &coord, std::array< double, D > &real_coord)
Given a path info object, get the coordinate from a given index into the patch.
Definition: DomainTools.h:134
ThunderEgg::DomainTools
Various tools for filling in values in a domain.
Definition: DomainTools.h:36
ThunderEgg::DomainTools::SetValues
static void SetValues(const Domain< 3 > &domain, Vector< 3 > &vec, int component_index, std::function< double(double, double, double)> func)
Set the values for a vector with the given function.
Definition: DomainTools.h:246
ThunderEgg::DomainTools::SetValues
static void SetValues(const Domain< D > &domain, Vector< D > &vec, int component_index, std::function< double(const std::array< double,(int) D > &)> func)
Set the values for a vector with the given function.
Definition: DomainTools.h:220
ThunderEgg::PatchInfo::spacings
std::array< double, D > spacings
The cell spacings in each direction.
Definition: PatchInfo.h:125
ThunderEgg::PatchInfo::starts
std::array< double, D > starts
The lower-left-bottom index of the patch.
Definition: PatchInfo.h:121
ThunderEgg::PatchInfo::num_ghost_cells
int num_ghost_cells
Number of ghost cells on each side of the patch.
Definition: PatchInfo.h:103
ThunderEgg::DomainTools::SetValues
static void SetValues(const Domain< 1 > &domain, Vector< 1 > &vec, int component_index, std::function< double(double)> func)
Set the values for a vector with the given function.
Definition: DomainTools.h:310
ThunderEgg::Vector
Vector class for use in thunderegg.
Definition: Vector.h:42
ThunderEgg::Domain::getPatchInfoVector
const std::vector< PatchInfo< D > > & getPatchInfoVector() const
Get a vector of PatchInfo pointers where index in the vector corresponds to the patch's local index.
Definition: Domain.h:259
ThunderEgg::Face
Enum-style class for the faces of an n-dimensional cube.
Definition: Face.h:41
ThunderEgg::DomainTools::SetValuesWithGhost
static void SetValuesWithGhost(const Domain< 1 > &domain, Vector< 1 > &vec, int component_index, std::function< double(double)> func)
Set the values (including ghost values) for a vector with the given function.
Definition: DomainTools.h:449
ThunderEgg::DomainTools::GetRealCoordBound
static void GetRealCoordBound(const PatchInfo< D > &pinfo, const std::array< int, D - 1 > &coord, Side< D > s, std::array< double, D > &real_coord)
Given a path info object and a side of the patch, get the coordinate from a given index into the inte...
Definition: DomainTools.h:178
ThunderEgg::Vector::getNumLocalPatches
int getNumLocalPatches() const
Get the number of local patches.
Definition: Vector.h:316
ThunderEgg::DomainTools::Integrate
static double Integrate(const Domain< D > &domain, const Vector< D > &u)
Integrate a vector over the domain.
Definition: DomainTools.h:513
ThunderEgg::DomainTools::SetValuesWithGhost
static void SetValuesWithGhost(const Domain< 3 > &domain, Vector< 3 > &vec, std::function< double(double, double, double)> func, Args... args)
Set the values (including ghost values) for a vector with the given functions.
Definition: DomainTools.h:499
ThunderEgg::View::getGhostEnd
const std::array< int, D > & getGhostEnd() const
Get the coordinate of the last ghost cell element.
Definition: View.h:219
ThunderEgg::View::getStart
const std::array< int, D > & getStart() const
Get the coordinate of the first element.
Definition: View.h:207
ThunderEgg::RuntimeError
ThunderEgg runtime exception.
Definition: RuntimeError.h:36
ThunderEgg::DomainTools::SetValues
static void SetValues(const Domain< 2 > &domain, Vector< 2 > &vec, int component_index, std::function< double(double, double)> func)
Set the values for a vector with the given function.
Definition: DomainTools.h:280
ThunderEgg::DomainTools::SetValuesWithGhost
static void SetValuesWithGhost(const Domain< D > &domain, Vector< D > &vec, std::function< double(const std::array< double, D > &)> func, Args... args)
Set the values (including ghost values) for a vector with the given functions.
Definition: DomainTools.h:481