sandbox/Antoonvh/topdown-iterator.h
We write grid-cell iterators the go from top to bottom, rather than the default bottom to top. This could be useful in exotic scenarios.
Normal N-order for trees:
2 4 | | 1 3
new mirrored-N order:
1 3 | / | 2 4
#include "foreach_cell_usd.h"
Now we can create usd versions of foreach_level(int l)
and foreach_level_or_leaf()
and foreach_usd()
, which concludes the excersize.
Grid * grid2 = NULL;
#define tree2 ((Tree *) grid2)
@def foreach_level2(l) {
if (l <= depth()) {
update_cache_f2();
CacheLevel _active = tree2->active[l];
foreach_cache_level (_active,l)
@
@define end_foreach_level2() end_foreach_cache_level(); }}
@define foreach_coarse_level2(l) foreach_level2(l) if (!is_leaf(cell)) {
@define end_foreach_coarse_level2() } end_foreach_level2()
@def foreach_level_or_leaf2(l) {
for (int _l1 = l; _l1 >= 0; _l1--)
foreach_level2(_l1)
if (_l1 == l || is_leaf (cell)) {
@
@define end_foreach_level_or_leaf2() } end_foreach_level2(); }
@define foreach_usd() update_cache_f2(); foreach_cache(tree2->leaves)
@define end_foreach_usd() end_foreach_cache()
A function that Updates the caches of grid2
.
trace
static void update_cache_f2 (void)
{
grid2 = grid;
Tree * q = tree2;
check_periodic (q);
foreach_cache (q->vertices)
if (level <= depth() && allocated(0))
cell.flags &= ~vertex;
/* empty caches */
q->leaves.n = q->faces.n = q->vertices.n = 0;
for (int l = 0; l <= depth(); l++)
q->active[l].n = q->prolongation[l].n =
q->boundary[l].n = q->restriction[l].n = 0;
#if FBOUNDARY
const unsigned short fboundary = 1 << user;
foreach_cell2() {
#else
foreach_cell_all2() {
#endif
if (is_local(cell) && is_active(cell)) {
// active cells
// assert (is_active(cell));
cache_level_append (&q->active[level], point);
}
#if !FBOUNDARY
if (is_boundary(cell)) {
// boundary conditions
bool has_neighbors = false;
foreach_neighbor (BGHOSTS)
if (allocated(0) && !is_boundary(cell))
has_neighbors = true, break;
if (has_neighbors)
cache_level_append (&q->boundary[level], point);
// restriction for masked cells
if (level > 0 && is_local(aparent(0)))
cache_level_append (&q->restriction[level], point);
}
#else
// boundaries
if (!is_boundary(cell)) {
// look in a 5x5 neighborhood for boundary cells
foreach_neighbor (BGHOSTS)
if (allocated(0) && is_boundary(cell) && !(cell.flags & fboundary)) {
cache_level_append (&q->boundary[level], point);
cell.flags |= fboundary;
}
}
// restriction for masked cells
else if (level > 0 && is_local(aparent(0)))
cache_level_append (&q->restriction[level], point);
#endif
if (is_leaf (cell)) {
if (is_local(cell)) {
cache_append (&q->leaves, point, 0);
// faces
unsigned short flags = 0;
foreach_dimension()
if (is_boundary(neighbor(-1)) || is_prolongation(neighbor(-1)) ||
is_leaf(neighbor(-1)))
flags |= face_x;
if (flags)
cache_append (&q->faces, point, flags);
foreach_dimension()
if (is_boundary(neighbor(1)) || is_prolongation(neighbor(1)) ||
(!is_local(neighbor(1)) && is_leaf(neighbor(1))))
cache_append (&q->faces, neighborp(1), face_x);
// vertices
for (int i = 0; i <= 1; i++)
#if dimension >= 2
for (int j = 0; j <= 1; j++)
#endif
#if dimension >= 3
for (int k = 0; k <= 1; k++)
#endif
if (!is_vertex(neighbor(i,j,k))) {
cache_append (&q->vertices, neighborp(i,j,k), 0);
neighbor(i,j,k).flags |= vertex;
}
// halo prolongation
if (cell.neighbors > 0)
cache_level_append (&q->prolongation[level], point);
}
else if (!is_boundary(cell) || is_local(aparent(0))) { // non-local
// faces
unsigned short flags = 0;
foreach_dimension()
if (allocated(-1) &&
is_local(neighbor(-1)) && is_prolongation(neighbor(-1)))
flags |= face_x;
if (flags)
cache_append_face (point, flags);
foreach_dimension()
if (allocated(1) && is_local(neighbor(1)) &&
is_prolongation(neighbor(1)))
cache_append_face (neighborp(1), face_x);
}
#if FBOUNDARY // fixme: this should always be included
continue;
#endif
}
}
/* optimize caches */
cache_shrink (&q->leaves);
cache_shrink (&q->faces);
cache_shrink (&q->vertices);
for (int l = 0; l <= depth(); l++) {
cache_level_shrink (&q->active[l]);
cache_level_shrink (&q->prolongation[l]);
cache_level_shrink (&q->boundary[l]);
cache_level_shrink (&q->restriction[l]);
}
q->dirty = false;
#if FBOUNDARY
for (int l = depth(); l >= 0; l--)
foreach_boundary_level (l)
cell.flags &= ~fboundary;
#endif
// mesh size
grid->n = q->leaves.n;
// for MPI the reduction operation over all processes is done by balance()
@if !_MPI
grid->tn = grid->n;
grid->maxdepth = grid->depth;
@endif
}