sandbox/Antoonvh/my_vertex.h
My Vertex methods
For finite-difference applications on trees, the requirement of an exact restriction operator warrants the use of vertex-points instead of cell centers. Furtheromore, for 5-point-stencil schemes, not all “Basilisk vertices” can be vertices.
User-interface vertex iterator
Iterate leafs and vertices on levels
@def foreach_vert()
update_cache();
foreach_cache(tree->leaves) {
x -= Delta/2.;
#if dimension >= 2
y -= Delta/2.;
#endif
#if dimension >= 3
z -= Delta/2.;
#endif
@
@define end_foreach_vert() } end_foreach_cache()
@def foreach_vert_level(l) {
if (l <= depth()) {
update_cache();
CacheLevel _active = tree->active[l];
foreach_cache_level (_active,l) {
x -= Delta/2.;
#if dimension >= 2
y -= Delta/2.;
#endif
#if dimension >= 3
z -= Delta/2.;
#endif
@
@define end_foreach_vert_level() } end_foreach_cache_level(); }}
Vertex restriction
Restrict only the local point
instead of 2^{\mathtt{dimension}}
static inline void restriction_vert (Point point, scalar s) {
s[] = fine(s,0,0,0);
}
Or a coarse estimate (not in 3D);
static inline void restriction_coarsen_vert (Point point, scalar s) {
#if (dimension == 1)
s[] = (fine(s,1,0,0) + 2*fine(s,0,0,0) + fine(s,-1,0,0))/4.;
#elif (dimension == 2)
s[] = (fine(s,1,0,0) + 2*fine(s,0,0,0) + fine(s,-1,0,0) +
fine(s,0,1,0) + fine(s,0,-1,0))/6.;
#endif
}
High-Level Boundary treatment
2nd-order
static inline void refine_vert (Point point, scalar s) {
// Injection
fine (s,0,0,0) = s[];
// Vertices with two nearest coarse neighbors
fine (s,1,0,0) = (s[] + s[1])/2.;
#if dimension > 1
fine (s,0,1,0) = (s[] + s[0,1])/2.;
#endif
#if dimension > 2
fine (s,0,0,1) = (s[] + s[0,0,1])/2.;
#endif
// Vertices with four nearest coarse neighbors
#if dimension > 1
fine(s,1,1,0) = (s[0] + s[1] + s[0,1] + s[1,1])/4.;
#endif
#if dimension > 2 // dimension == 3
fine(s,1,0,1) = (s[] + s[1] + s[0,0,1] + s[1,0,1])/4.;
fine(s,0,1,1) = (s[] + s[0,1] + s[0,1] + s[0,1,1])/4.;
// In 3D, there is a vertex with 8 nearest coarse neighbors
fine(s,1,1,1) = (s[] + s[1,1,1] +
s[1] + s[0,1] + s[0,0,1] +
s[1,1] + s[0,1,1] + s[1,0,1])/8.;
#endif
}
static inline void prolongation_vert (Point point, scalar s) {
// do not inject for error estimation
refine_vert(point, s);
//if (!is_leaf(cell))
fine(s,0,0,0) = (s[-1] + s[1]
#if dimension > 1
+ s[0,1] + s[0,-1]
#endif
)/(2*dimension);
}
4th-order
static inline void refine_vert4 (Point point, scalar s) {
// Injection
fine (s,0,0,0) = s[];
// Vertices with two nearest coarse neighbors
fine (s,1,0,0) = (9*(s[0] + s[1]) - (s[-1] + s[2]))/16.;
#if dimension > 1
fine (s,0,1,0) = (9*(s[0] + s[0,1]) - (s[0,-1] + s[0,2]))/16.;
#endif
#if dimension > 2
fine (s,0,0,1) = (9*(s[0] + s[0,0,1]) - (s[0,0,-1] + s[0,0,2]))/16.;
#endif
// Vertices with four nearest coarse neighbors
#if dimension > 1
fine(s,1,1,0) = (9.*(s[0] + s[1,1] + s[0,1] + s[1,0]) -
(s[-1,-1] + s[2,2] + s[-1,2] + s[2,-1]))/32.;
#endif
// not 3D yet...
#if dimension > 2 // dimension == 3
fine(s,1,0,1) = (s[] + s[1] + s[0,0,1] + s[1,0,1])/4.;
fine(s,0,1,1) = (s[] + s[0,1] + s[0,1] + s[0,1,1])/4.;
// In 3D, there is a vertex with 8 nearest coarse neighbors
fine(s,1,1,1) = (s[] + s[1,1,1] +
s[1] + s[0,1] + s[0,0,1] +
s[1,1] + s[0,1,1] + s[1,0,1])/8.;
#endif
}
static inline void prolongation_vert4 (Point point, scalar s) {
refine_vert4(point, s);
if (!is_leaf(cell)) {
fine(s,0,0,0) = (4*(fine(s,-1,0,0) + fine(s,1,0,0)) -
(fine(s,-2,0,0) + fine(s,2,0,0))
#if dimension > 1
+4* (fine(s,0,-1,0) + fine(s,0,1,0)) -
(fine(s,0,-2,0) + fine(s,0,2,0))
#endif
)/(dimension*6);
}
}
foreach_dimension()
double interp4_x (Point point, scalar s) {
return (4*(s[-1] + s[1])- (s[-2] + s[2]))/6.;
}
Biassed 5th order
double interp5_x (Point point, scalar s) {
return ((3./128)*coarse(s,-2,0,0) - (5./32)*coarse(s,-1,0,0) +
(45./64.)*coarse(s,0,0,0) + (15./32.)*coarse(s,1,0,0) -
(5./128.)*coarse(s,2,0,0));
}
double interp5_y (Point point, scalar s) {
return ((3./128)*coarse(s,0,-2,0) - (5./32)*coarse(s,0,-1,0) +
(45./64.)*coarse(s,0,0,0) + (15./32.)*coarse(s,0,1,0) -
(5./128.)*coarse(s,0,2,0));
}
static inline void refine_vert5 (Point point, scalar s) {
// Injection
fine (s,0,0,0) = s[];
// Vertices with two nearest coarse neighbors
fine (s,1,0,0) = ((3./128)*s[-2] - (5./32)*s[-1] + (45./64.)*s[]
+ (15./32.)*s[1] - (5./128.)*s[2]);
#if dimension > 1
fine (s,0,1,0) = ((3./128)*s[0,-2] - (5./32)*s[0,-1] + (45./64.)*s[]
+ (15./32.)*s[0,1] - (5./128.)*s[0,2]);
#endif
#if dimension > 2
fine (s,0,0,1) = ((3./128)*s[0,0,-2] - (5./32)*s[0,0,-1] + (45./64.)*s[]
+ (15./32.)*s[0,0,1] - (5./128.)*s[0,0,2]);
#endif
// Vertices with four nearest coarse neighbors
#if dimension > 1
fine (s,1,1,0) = ((3./128)*s[-2,-2] - (5./32)*s[-1,-1] + (45./64.)*s[] +
(15./32.)*s[1,1] - (5./128.)*s[2,2]);
#endif
#if dimension > 2 // dimension == 3
fine(s,1,0,1) = ((3./128)*s[-2,0,-2] - (5./32)*s[-1,0,-1] + (45./64.)*s[] +
(15./32.)*s[1,0,1] - (5./128.)*s[2,0,2]);
fine(s,0,1,1) = ((3./128)*s[0,-2,-2] - (5./32)*s[0,-1,-1] + (45./64.)*s[] +
(15./32.)*s[0,1,1] - (5./128.)*s[0,2,2]);
// In 3D, there is a vertex with 8 nearest coarse neighbors
fine(s,1,1,1) = ((3./128)*s[-2,-2,-2] - (5./32)*s[-1,-1,-1] + (45./64.)*s[] +
(15./32.)*s[1,1,1] - (5./128.)*s[2,2,2]);
#endif
}
To do
- Other than periodic box boundaries