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
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
| #include "output_xdmf.h"
#define r2 (sq(x) + sq(y))
int main()
{
L0 = 1. [1];
origin(-L0 / 2, -L0 / 2, -L0 / 2);
N = 1 << (MAXLEVEL - 3);
init_grid(N);
#if TREE
double outer_radius = 0.25 [1];
double inner_radius = 0.1 [1];
refine(((r2 < sq(outer_radius)) && (r2 > sq(inner_radius))) && level < MAXLEVEL);
#endif
scalar s[], c[], f[];
vector u[];
coord k = {1. [-1], 1. [-1], (4. * pi)[-1]};
double yref1 = 0.1 [1];
double yref2 = -0.25 [1];
double thck = 0.02 [1];
foreach ()
{
s[] = sin(k.x * x) * cos(k.y * y) * cos(k.z * z) + 0.1 * noise();
f[] = (y < yref1);
c[] = 0.5 * (1. + tanh((yref2 - y) / thck));
foreach_dimension()
u.x[] = sq(x) + sq(y);
}
boundary({s, c, u});
int num_points = 0, num_points_loc=0;
foreach_vertex(reduction(+:num_points)) {
num_points++;
}
foreach_vertex(serial, noauto) {
num_points_loc++;
}
int num_cells = 0, num_cells_loc=0;
foreach(reduction(+:num_cells)) {
num_cells++;
}
foreach(serial, noauto) {
num_cells_loc++;
}
fprintf(stdout, "%d %d %d %d %d \n", pid(), num_points, num_points_loc, num_cells, num_cells_loc);
output_xmf({s, f, c}, {u}, "domain");
#if dimension > 2
output_xmf_slice({s, f, c}, {u}, "slice_x", (coord){1, 0, 0}, 0);
output_xmf_slice({s, f, c}, {u}, "slice_y", (coord){0, 1, 0}, 0);
output_xmf_slice({s, f, c}, {u}, "slice_z", (coord){0, 0, 1}, 0);
#endif
}
|
