/*  prog1a.c -- Solve Poisson Equation
    Written Feb 25, 2003 by Eric Olson for Mathematics 484

    For information about the five-point formula for solving the Poisson
    equation see Chapter 7 in A First Course in the Numerical Analysis
    of Differential Equations by Iserles.  For information about the
    Gauss-Seidel iteration method used to solve the resulting system of
    finite differences see Chapter 10 in the same reference.  */

#include <stdio.h>
#include <math.h>

#define N 10
const double h=8.0/N;
const double eps=0.001;
const double etol=0.001/N/N/N;
double u[N][N],hhf[N][N];
int cmax=0;
struct intpt {
    int i,j,n;
} ipts[2*N];

/*  Determine if a point (x,y) is an interior point in the domain.
    Returns 1 if an interior point 0 otherwise */
int isint(double x, double y){
    if( (x>0.0+eps) &&
        (y>0.0+eps) &&
        (x+y<7.0-eps) ) return 1;
    else return 0;
}

/*  Specify the Dirichlet condition at (x,y) on the boundary */
double psi(double x, double y){
    return x-2.0*y;
}

/*  Specify the forcing term at (x,y) in the domain */
double f(double x, double y){
    return 0.0;
}

void printu(){
    int i,j;
    /*  Print the matrix u */
    for(i=0;i<N;i++){
        for(j=0;j<N;j++){
            printf("%5.1f ",u[i][j]);
        }
        printf("\n");
    }
}

void doinit(){
    int i,j,c;

    /*  Locate the domain */
    for(i=0,c=0;i<N;i++){
        ipts[c].n=0;
            for(j=0;j<N;j++){
            if(isint(i*h,j*h)){
                u[i][j]=1.0;
                if(!ipts[c].n){
                    ipts[c].i=i;
                    ipts[c].j=j;
                }
                ipts[c].n++;
            } else {
                u[i][j]=0.0;
                if(ipts[c].n){
                    c++;
                    ipts[c].n=0;
                }
            }
        }
    }
    cmax=c;

    /*  Initialize the exterior points */
    for(c=0;c<cmax;c++){
        int jmax=ipts[c].j+ipts[c].n;
        i=ipts[c].i;
        for(j=ipts[c].j;j<jmax;j++){
            if(!u[i+1][j]) u[i+1][j]=psi((i+1)*h,j*h);
            if(!u[i-1][j]) u[i-1][j]=psi((i-1)*h,j*h);
            if(!u[i][j+1]) u[i][j+1]=psi(i*h,(j+1)*h);
            if(!u[i][j-1]) u[i][j-1]=psi(i*h,(j-1)*h);
        }
    }

    /*  Initialize the interior and forcing points */
    for(c=0;c<cmax;c++){
        int jmax=ipts[c].j+ipts[c].n;
        i=ipts[c].i;
        for(j=ipts[c].j;j<jmax;j++){
            u[i][j]=0.0;
            hhf[i][j]=h*h*f(i*h,j*h);
        }
    }
}

void gauss_seidel(){
    /*  Solve the problem */
    int n;
    for(n=0;n<100000;n++){
        int i,j,c;
        double uerr=0.0;
        for(c=0;c<cmax;c++){
            int jmax=ipts[c].j+ipts[c].n;
            i=ipts[c].i;
            for(j=ipts[c].j;j<jmax;j++){
                double unew=(u[i][j-1]+u[i][j+1]+u[i+1][j]
                    +u[i-1][j]-hhf[i][j])/4.0;
                double udif=unew-u[i][j];
                u[i][j]=unew;
                uerr+=udif*udif;
            }
        }
        if(uerr<=etol*etol) break;
    }
    printf("Convergence after %d iterations\n",n);
}

main(){
    doinit();
    gauss_seidel();
    printu();
    exit(0);
}