src/test/boundary_vertex3D.c
Checks boundary conditions for vertex fields in parallel (in 3D).
#include "grid/octree.h"
#define shape(x,y,z) (sq(x) + sq(y) + sq(z))
int main()
{
size (16);
init_grid (4);
origin (-L0/2., -L0/2., -L0/2.);
refine (sq(x) + sq(y) + sq(z) < sq(L0/4.) && level < 3);
scalar s[];
foreach_dimension() {
s[left] = dirichlet(shape(x,y,z));
s[right] = dirichlet(shape(x,y,z));
}
foreach()
s[] = shape(x,y,z);
vertex scalar phi[];
foreach_vertex()
phi[] = (s[] + s[-1] + s[0,-1] + s[-1,-1] +
s[0,0,-1] + s[-1,0,-1] + s[0,-1,-1] + s[-1,-1,-1])/8.;
foreach() {
fprintf (qerr, "%g %g %g %g\n", x - Delta/3., y - Delta/3., z - Delta/3.,
phi[]);
fprintf (qerr, "%g %g %g %g\n", x - Delta/3., y + Delta/3., z - Delta/3.,
phi[0,1]);
fprintf (qerr, "%g %g %g %g\n", x + Delta/3., y - Delta/3., z - Delta/3.,
phi[1,0]);
fprintf (qerr, "%g %g %g %g\n", x + Delta/3., y + Delta/3., z - Delta/3.,
phi[1,1]);
fprintf (qerr, "%g %g %g %g\n", x - Delta/3., y - Delta/3., z + Delta/3.,
phi[0,0,1]);
fprintf (qerr, "%g %g %g %g\n", x - Delta/3., y + Delta/3., z + Delta/3.,
phi[0,1,1]);
fprintf (qerr, "%g %g %g %g\n", x + Delta/3., y - Delta/3., z + Delta/3.,
phi[1,0,1]);
fprintf (qerr, "%g %g %g %g\n", x + Delta/3., y + Delta/3., z + Delta/3.,
phi[1,1,1]);
}
#if _MPI
debug_mpi (NULL);
#else
output_cells();
#endif
}