sandbox/bderembl/filter.c
Filter function
#include "grid/multigrid.h"
#include "run.h"
The function below will replace the default one.
void waveletb (scalar s, scalar w)
{
restriction ({s});
for (int l = depth() - 1; l >= 0; l--) {
foreach_coarse_level (l) {
foreach_child()
w[] = s[];
s.prolongation (point, s);
foreach_child() {
double sp = s[];
s[] = w[];
/* difference between fine value and its prolongation */
w[] -= sp;
}
}
boundary_level ({w}, l + 1);
}
/* root cell */
foreach_level(0)
w[] = s[];
boundary_level ({w}, 0);
}
Inverse wavelet transform
void inverse_waveletb (scalar s, scalar w)
{
foreach_level(0)
s[] = w[];
boundary_level ({s}, 0);
for (int l = 0; l <= depth() - 1; l++) {
foreach_coarse_level (l) {
s.prolongation (point, s);
foreach_child()
s[] += w[];
}
boundary_level ({s}, l + 1);
}
}
Main loop
int main() {
N = 1 << 8;
init_grid (N);
scalar s[], w[], s2[];
s[top] = 0;
s[bottom] = 0;
s[right] = 0;
s[left] = 0;
w[top] = 0;
w[bottom] = 0;
w[right] = 0;
w[left] = 0;
s2[top] = 0;
s2[bottom] = 0;
s2[right] = 0;
s2[left] = 0;
foreach()
s[] = sin(2*pi*x)*sin(10*pi*y);
boundary({s});
Wavelet transform
waveletb (s, w);
inverse_waveletb (s2, w);
Check that we recover the original signal (almost) exactly.
foreach()
assert (fabs(s[] - s2[]) < 1e-10);
char name[80];
sprintf (name,"out.dat");
FILE * fp = fopen (name, "w");
output_field ({s, s2}, fp, linear=1);
fclose(fp);
}
We plot the filtered field and the filter
import numpy as np
import matplotlib.pyplot as plt
file1 = 'out.dat'
data = np.loadtxt(file1)
n1, n2 = data.shape
N = int(np.sqrt(n1))
x = data[:,0].reshape((N,N))
y = data[:,1].reshape((N,N))
p0 = data[:,2].reshape((N,N))
p1 = data[:,3].reshape((N,N))
plt.figure()
plt.contourf(x,y,p0)
plt.savefig('field.png')
plt.figure()
plt.contourf(x,y,p1)
plt.savefig('filter.png')