#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<malloc.h>
#include <iostream>
#include <string>
using namespace std;
#define e_eps 1.0e-8
#define PI 3.14159265358979323846264338327950288
 
void Set_kyokai(double **U,int Nx,int Ny){
        int x,y;
        double left, top, right, bottom;
 
        cout<<"境界条件\n";
        cout<<"left :";
        cin>>left;
        cout<<"right :";
        cin>>right;
        cout<<"bottom :";
        cin>>bottom;
 
        double Hx = 1.0/ (double)Nx;
        double Hy = 1.0/(double)Ny;
 
        /* 領域の底の境界 */
        for (x = 0; x <= Nx; x++){
                U[x][0] = bottom;
                printf("%lf\n",U[x][0]);
        }
        cout<<"\n";
 
        /* 領域の上の境界 */
        for (x = 0; x <= Nx; x++){
                top = sin(PI * Hx *x);
                U[x][Ny] = top;
                printf("%lf\n",U[x][Ny]);
        }
 
        /* 左の壁の境界 */
        for (y = 0; y <= Ny; y++){
                U[0][y] = left;
        }
 
        /* 右の壁の境界 */
        for (y = 0; y <= Ny; y++){
                U[Nx][y] = right;
        }
}
 
int main(void){
        int i,j,k;
        double uu = 0.0, d = 0.0, sum = 0.0;
        double errMAX;
 
        int Nx,Ny;
        Nx = Ny = 0;
 
        //計算領域の分割数
        cout<<"分割数の設定\n";
        cout<<"x軸 N :";
        cin>>Nx;
        cout<<"y軸 N :";
        cin>>Ny;
 
        //緩和係数
        double ω = (2.0)/(1.0 + sin(PI/Nx));
        
        double **D = new double* [Nx+1];
        for (i = 0; i < Ny +1; i++) {
                D[i] = new double[Ny+1];
        }
 
        Set_kyokai(D,Nx,Ny);
 
        for(k=1;k<=1000;k++){
                errMAX = 0.0;
                sum = 0.0;
                for (i = 1; i <= Nx -1; i++){
                        for (j = 1; j <= Ny -1; j++){
                                /* SORによる解法 */                                
                                uu = D[i+1][j]+D[i-1][j]+D[i][j+1]+D[i][j-1];
                                d = uu * ω/4.0 + (1.0 - ω)*D[i][j];
 
                                sum += fabs(d);
 
                                errMAX += fabs(d - D[i][j]);
                                D[i][j] = d;
                        }
                }
                if(errMAX < e_eps * sum){
                        break;
                }
        }
 
        cout<<"結果↓\n";
        for(i=0;i<Nx +1;i++){
                for(j=0;j<Ny +1;j++){
                        printf("U[%d][%d] =%9f\n",i,j,D[i][j]);
                }
                cout<<"\n";
        }
        printf("\n");
        
        for (i = 0; i < Ny +1; i++) {
                delete[] D[i];
        }
        delete[] D;
 
        return 0;
}