// #include "grid/multigrid.h"
#include "utils.h"
#include "poisson.h"

double solution(double x, double y) { return (sin(pi * x) + cos(pi * y)); }
scalar U[], f[], ee[];
// int nrelax = 10000;
int depth = 8;
// void * data;
struct Poisson pp;
const scalar lambda[] = 0;

U[right] = dirichlet(solution(x, y));
U[left] = dirichlet(solution(x, y));
U[top] = dirichlet(solution(x, y));
U[bottom] = dirichlet(solution(x, y));

int main()
{
    pp.a = U;
    pp.b = f;
    pp.alpha = unityf;
    pp.lambda = lambda;
    L0 = 2 * pi;
    X0 = -L0 / 2.;
    Y0 = -L0 / 2.;

    init_grid(1 << depth);
    foreach ()
    {
        U[] = 0;
        f[] = -sq(pi) * sin(pi * x) - sq(pi) * cos(pi * y);
    }
    double Tolerance = 2;
    double sumerror = 0.;
    foreach ()
    {
        sumerror += fabs(solution(x, y));
    }
    printf("%g\n", sumerror);
    int j = 0;
    while (sumerror > Tolerance)
    {
        double e = 0;
        relax({U}, {f}, depth, &pp);
        foreach ()
        {
            e += fabs(U[] - solution(x, y));
        }
        if (e < sumerror)
            sumerror = e;
        printf("%d %g %g\n", j, sumerror, e);
        j++;
    }

    foreach ()
    {
        ee[] = U[] - solution(x, y);
    }
    FILE *fp1 = fopen("outputfield_U", "w");
    output_field({U}, fp1);
    FILE *fp2 = fopen("outputfield_e", "w");
    output_field({ee}, fp2);
}