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
| #include "output_vtu.h"
#define r2 (sq(x) + sq(y))
int main(){
L0 = 1. [1];
origin(-L0 / 2, -L0 / 2, -L0 / 2);
N = 1<<(MAXLEVEL-4);
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});
output_vtu({s,f,c}, {u}, "domain");
#if dimension > 2
output_slice_vtu({s,f,c}, {u}, "slice_x", (coord){1,0,0}, 0);
output_slice_vtu({s,f,c}, {u}, "slice_z");
#endif
}
|
