1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
| /* similar to circle.c but in 3D (see also poisson3D.c) */
#include "grid/octree.h"
#include "utils.h"
#include "poisson.h"
double solution (double x, double y, double z)
{
return sin(3.*pi*x)*sin(3.*pi*y)*sin(3.*pi*z);
}
int main()
{
size (1.[0]);
origin (-0.5, -0.5, -0.5);
int depth = 6;
init_grid(1);
refine (level < depth - 2 || level <= depth*(1. - sqrt(x*x + y*y + z*z)));
scalar a[], b[];
foreach_dimension() {
a[right] = dirichlet (solution(x, y, z));
a[left] = dirichlet (solution(x, y, z));
}
foreach() {
a[] = 0.;
b[] = -27.*pi*pi*sin(3.*pi*x)*sin(3.*pi*y)*sin(3.*pi*z);
}
poisson (a, b);
double max = 0;
foreach(reduction(max:max)) {
double e = a[] - solution(x, y, z);
if (fabs(e) > max) max = fabs(e);
printf ("%g %g %g %g %g %g\n", x, y, z, a[], b[], e);
}
fprintf (stderr, "max error %g\n", max);
}
|
