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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
| struct Adapt2 {
scalar * slist; // list of scalars
double * max; // tolerance for each scalar
int * maxlevel; // maximum level of refinement for each scalar
int minlevel; // minimum level of refinement (default 1)
scalar * list; // list of fields to update (default all)
};
trace
astats adapt_wavelet2 (struct Adapt2 p)
{
if (p.list == NULL)
p.list = all;
if (is_constant(cm))
restriction (p.slist);
else {
scalar * listr = list_concat ({cm}, p.slist);
restriction (listr);
free (listr);
}
astats st = {0, 0};
scalar * listc = NULL;
for (scalar s in p.list)
if (!is_constant(s) && s.restriction != no_restriction)
listc = list_add (listc, s);
// refinement
if (p.minlevel < 1)
p.minlevel = 1;
tree->refined.n = 0;
static const int refined = 1 << user, too_fine = 1 << (user + 1);
foreach_cell() {
if (is_active(cell)) {
static const int too_coarse = 1 << (user + 2);
if (is_leaf (cell)) {
if (cell.flags & too_coarse) {
cell.flags &= ~too_coarse;
refine_cell (point, listc, refined, &tree->refined);
st.nf++;
}
continue;
}
else { // !is_leaf (cell)
if (cell.flags & refined) {
// cell has already been refined, skip its children
cell.flags &= ~too_coarse;
continue;
}
// check whether the cell or any of its children is local
bool local = is_local(cell);
if (!local)
foreach_child()
if (is_local(cell))
local = true, break;
if (local) {
int i = 0;
static const int just_fine = 1 << (user + 3);
for (scalar s in p.slist) {
double max = p.max[i], sc[1 << dimension];
int mlev = p.maxlevel[i++];
int c = 0;
foreach_child()
sc[c++] = s[];
s.prolongation (point, s);
c = 0;
foreach_child() {
double e = fabs(sc[c] - s[]);
if (e > max && level < mlev) { //tag cells that are too coarse
cell.flags &= ~too_fine;
cell.flags |= too_coarse;
}
else if ((e <= max/1.5 || level > (mlev-1)) &&
!(cell.flags & (too_coarse|just_fine))) { //tag cells that are too fine
if (level >= p.minlevel)
cell.flags |= too_fine;
}
else if (!(cell.flags & too_coarse)) { // Tag cells that are just fine
cell.flags &= ~too_fine;
cell.flags |= just_fine;
}
s[] = sc[c++];
}
}
foreach_child() {
cell.flags &= ~just_fine;
if (!is_leaf(cell)) {
cell.flags &= ~too_coarse;
if (level >= 13)
cell.flags |= too_fine;
}
else if (!is_active(cell))
cell.flags &= ~too_coarse;
}
}
}
}
else // inactive cell
continue;
}
mpi_boundary_refine (listc);
// coarsening
// the loop below is only necessary to ensure symmetry of 2:1 constraint
for (int l = depth(); l >= 0; l--) {
foreach_cell()
if (!is_boundary(cell)) {
if (level == l) {
if (!is_leaf(cell)) {
if (cell.flags & refined)
// cell was refined previously, unset the flag
cell.flags &= ~(refined|too_fine);
else if (cell.flags & too_fine) {
if (is_local(cell) && coarsen_cell (point, listc))
st.nc++;
cell.flags &= ~too_fine; // do not coarsen parent
}
}
if (cell.flags & too_fine)
cell.flags &= ~too_fine;
else if (level > 0 && (aparent(0).flags & too_fine))
aparent(0).flags &= ~too_fine;
continue;
}
else if (is_leaf(cell))
continue;
}
mpi_boundary_coarsen (l, too_fine);
}
free (listc);
mpi_all_reduce (st.nf, MPI_INT, MPI_SUM);
mpi_all_reduce (st.nc, MPI_INT, MPI_SUM);
if (st.nc || st.nf)
mpi_boundary_update (p.list);
return st;
}
|
