DTKQAPJB6JCJZYLQ4MOFA5QHDCCIQ3KG5ED6PVF5PGEH33YSJFBAC
ER6QNVNNYMJOPQ3UH3KONKGUIYLCPCR5AIAYD6BVKV43TDU6ZBUQC
YIQN7NJTGEVKW7JZHL6CTH6EPCIXCNBYNURIGXPYZAOUX3VAJQMAC
T3TZRPPAIA24I3YL3JFB4XEAYCWU3HJAJUCF7NNIFMP4I5X4SM5QC
FEMASUBNU32NSG4DNXZX54CGCA57PVRGYO46L3A6F2EJ4BCSJ3SAC
BJDGFYBMECTJG7BHLNHLSCUCBVYHAY6OGY37FIJP6JDGNDXQNQVAC
OHAARXLBTFOIV5FOTY7IDMOIAE275CKXCE6XQMWV3EC63YCSMIWAC
PG2P3JQPLWZRLCAJ7JY2B7G7AM5DHOE2CJUKIPWREI3NAUKZF3TAC
const auto kernel{[&](const int i, const int j,
const int k) CCTK_ATTRIBUTE_ALWAYS_INLINE {
const CCTK_REAL x = x0[0] + (lbnd[0] + i + CCTK_REAL(CI - 1) / 2) * dx[0];
const CCTK_REAL y = x0[1] + (lbnd[1] + j + CCTK_REAL(CJ - 1) / 2) * dx[1];
const CCTK_REAL z = x0[2] + (lbnd[2] + k + CCTK_REAL(CK - 1) / 2) * dx[2];
const int idx = i * di + j * dj + k * dk;
const vect<int, dim> I{{i, j, k}};
const PointDesc p{i, j, k, x, y, z, idx, dj, dk, I, inormal};
f(p);
}};
array<bool, dim> bforward;
for (int d = 0; d < dim; ++d)
bforward[d] = inormal[d] >= 0;
bool all_forward = true;
for (int d = 0; d < dim; ++d)
all_forward &= bforward[d];
if (all_forward) {
for (int k = imin[2]; k < imax[2]; ++k) {
for (int j = imin[1]; j < imax[1]; ++j) {
#pragma omp simd
for (int i = imin[0]; i < imax[0]; ++i) {
kernel(i, j, k);
}
}
}
for (int i = imin[0]; i < imax[0]; ++i) {
CCTK_REAL x = x0[0] + (lbnd[0] + i + CCTK_REAL(CI - 1) / 2) * dx[0];
CCTK_REAL y = x0[1] + (lbnd[1] + j + CCTK_REAL(CJ - 1) / 2) * dx[1];
CCTK_REAL z = x0[2] + (lbnd[2] + k + CCTK_REAL(CK - 1) / 2) * dx[2];
int idx = i * di + j * dj + k * dk;
vect<int, dim> I{{i, j, k}};
const PointDesc p{i, j, k, x, y, z, idx, dj, dk, I, inormal};
f(p);
for (int i0 = 0; i0 < imax[0] - imin[0]; ++i0) {
const int i = bforward[0] ? imin[0] + i0 : imax[0] - 1 - i0;
const int j = bforward[1] ? imin[1] + j0 : imax[1] - 1 - j0;
const int k = bforward[2] ? imin[2] + k0 : imax[2] - 1 - k0;
kernel(i, j, k);
}
for (int nk = -1; nk <= +1; ++nk) {
for (int nj = -1; nj <= +1; ++nj) {
for (int ni = -1; ni <= +1; ++ni) {
if ((ni != 0 && bbox[0 + (ni == -1 ? 0 : 1)]) ||
(nj != 0 && bbox[2 + (nj == -1 ? 0 : 1)]) ||
(nk != 0 && bbox[4 + (nk == -1 ? 0 : 1)])) {
for (int rank = dim - 1; rank >= 0; --rank) {
for (int nk = -1; nk <= +1; ++nk) {
for (int nj = -1; nj <= +1; ++nj) {
for (int ni = -1; ni <= +1; ++ni) {
if ((ni == 0) + (nj == 0) + (nk == 0) == rank) {
if ((ni != 0 && bbox[0 + (ni == -1 ? 0 : 1)]) ||
(nj != 0 && bbox[2 + (nj == -1 ? 0 : 1)]) ||
(nk != 0 && bbox[4 + (nk == -1 ? 0 : 1)])) {
array<int, dim> imin, imax;
for (int d = 0; d < dim; ++d) {
const int ghost_offset = nghostzones[d] - group_nghostzones[d];
const int begin_bnd = ghost_offset;
const int begin_int = nghostzones[d];
const int end_int = lsh[d] + offset[d] - nghostzones[d];
const int end_bnd = lsh[d] + offset[d] - ghost_offset;
switch (inormal[d]) {
case -1: // lower boundary
imin[d] = begin_bnd;
imax[d] = begin_int;
break;
case 0: // interior
imin[d] = begin_int;
imax[d] = end_int;
break;
case +1: // upper boundary
imin[d] = end_int;
imax[d] = end_bnd;
break;
default:
assert(0);
}
array<int, dim> imin, imax;
for (int d = 0; d < dim; ++d) {
const int ghost_offset =
nghostzones[d] - group_nghostzones[d];
const int begin_bnd = ghost_offset;
const int begin_int = nghostzones[d];
const int end_int = lsh[d] + offset[d] - nghostzones[d];
const int end_bnd = lsh[d] + offset[d] - ghost_offset;
switch (inormal[d]) {
case -1: // lower boundary
imin[d] = begin_bnd;
imax[d] = begin_int;
break;
case 0: // interior
imin[d] = begin_int;
imax[d] = end_int;
break;
case +1: // upper boundary
imin[d] = end_int;
imax[d] = end_bnd;
break;
default:
assert(0);
}
imin[d] = std::max(tmin[d], imin[d]);
imax[d] = std::min(tmax[d] + (tmax[d] >= lsh[d] ? offset[d] : 0),
imax[d]);
}
imin[d] = std::max(tmin[d], imin[d]);
imax[d] = std::min(
tmax[d] + (tmax[d] >= lsh[d] ? offset[d] : 0), imax[d]);
}
for (int nk = -1; nk <= +1; ++nk) {
for (int nj = -1; nj <= +1; ++nj) {
for (int ni = -1; ni <= +1; ++ni) {
if ((ni != 0 && !bbox[0 + (ni == -1 ? 0 : 1)]) ||
(nj != 0 && !bbox[2 + (nj == -1 ? 0 : 1)]) ||
(nk != 0 && !bbox[4 + (nk == -1 ? 0 : 1)])) {
for (int rank = dim - 1; rank >= 0; --rank) {
for (int nk = -1; nk <= +1; ++nk) {
for (int nj = -1; nj <= +1; ++nj) {
for (int ni = -1; ni <= +1; ++ni) {
if ((ni == 0) + (nj == 0) + (nk == 0) == rank) {
if ((ni != 0 && !bbox[0 + (ni == -1 ? 0 : 1)]) ||
(nj != 0 && !bbox[2 + (nj == -1 ? 0 : 1)]) ||
(nk != 0 && !bbox[4 + (nk == -1 ? 0 : 1)])) {
const array<int, dim> inormal{ni, nj, nk};
array<int, dim> imin, imax;
for (int d = 0; d < dim; ++d) {
const int ghost_offset =
nghostzones[d] - group_nghostzones[d];
const int begin_bnd = ghost_offset;
const int begin_int = nghostzones[d];
const int end_int = lsh[d] + offset[d] - nghostzones[d];
const int end_bnd = lsh[d] + offset[d] - ghost_offset;
switch (inormal[d]) {
case -1: // lower boundary
imin[d] = begin_bnd;
imax[d] = begin_int;
break;
case 0: // interior
imin[d] = begin_int;
imax[d] = end_int;
break;
case +1: // upper boundary
imin[d] = end_int;
imax[d] = end_bnd;
break;
default:
assert(0);
}
array<int, dim> imin, imax;
for (int d = 0; d < dim; ++d) {
const int ghost_offset = nghostzones[d] - group_nghostzones[d];
const int begin_bnd = ghost_offset;
const int begin_int = nghostzones[d];
const int end_int = lsh[d] + offset[d] - nghostzones[d];
const int end_bnd = lsh[d] + offset[d] - ghost_offset;
switch (inormal[d]) {
case -1: // lower boundary
imin[d] = begin_bnd;
imax[d] = begin_int;
break;
case 0: // interior
imin[d] = begin_int;
imax[d] = end_int;
break;
case +1: // upper boundary
imin[d] = end_int;
imax[d] = end_bnd;
break;
default:
assert(0);
loop_box<CI, CJ, CK>(f, imin, imax, inormal);
} // if rank
}
}
}
} // for rank
}
// Loop over all outer ghost points. This excludes ghost edges/corners on
// non-ghost faces. Loop over faces first, then edges, then corners.
template <int CI, int CJ, int CK, typename F>
void loop_ghosts(const array<int, dim> &group_nghostzones, const F &f) const {
constexpr array<int, dim> offset{!CI, !CJ, !CK};
imin[d] = std::max(tmin[d], imin[d]);
imax[d] = std::min(tmax[d] + (tmax[d] >= lsh[d] ? offset[d] : 0),
imax[d]);
}
for (int rank = dim - 1; rank >= 0; --rank) {
for (int nk = -1; nk <= +1; ++nk) {
for (int nj = -1; nj <= +1; ++nj) {
for (int ni = -1; ni <= +1; ++ni) {
if ((ni == 0) + (nj == 0) + (nk == 0) == rank) {
if ((ni == 0 || !bbox[0 + (ni == -1 ? 0 : 1)]) &&
(nj == 0 || !bbox[2 + (nj == -1 ? 0 : 1)]) &&
(nk == 0 || !bbox[4 + (nk == -1 ? 0 : 1)])) {
const array<int, dim> inormal{ni, nj, nk};
loop_box<CI, CJ, CK>(f, imin, imax, inormal);
array<int, dim> imin, imax;
for (int d = 0; d < dim; ++d) {
const int ghost_offset =
nghostzones[d] - group_nghostzones[d];
const int begin_bnd = ghost_offset;
const int begin_int = nghostzones[d];
const int end_int = lsh[d] + offset[d] - nghostzones[d];
const int end_bnd = lsh[d] + offset[d] - ghost_offset;
switch (inormal[d]) {
case -1: // lower boundary
imin[d] = begin_bnd;
imax[d] = begin_int;
break;
case 0: // interior
imin[d] = begin_int;
imax[d] = end_int;
break;
case +1: // upper boundary
imin[d] = end_int;
imax[d] = end_bnd;
break;
default:
assert(0);
}
imin[d] = std::max(tmin[d], imin[d]);
imax[d] = std::min(
tmax[d] + (tmax[d] >= lsh[d] ? offset[d] : 0), imax[d]);
}
loop_box<CI, CJ, CK>(f, imin, imax, inormal);
}
} // if rank