#include <iostream>
#include <cmath>
#include <algorithm>
#include <vector>

using namespace std;

#define A1 -2
#define A2 3
#define B1 -1
#define B2 4
#define M 20
#define N 20
#define SIZE (M + 1) * (N + 1)

const double h1 = (A2 - A1) * 1.0 / M;
const double h2 = (B2 - B1) * 1.0 / N;

vector<double> w;
vector<double> old_w;

// exact solution u(x, y)
double u(double x, double y) {
	return 2.0 / (1 + pow(x, 2) + pow(y, 2));
}

double k(double x, double y) {
	return 1 + pow(x + y, 2);
}

// potential function
double q(double x, double y) {
	return 1;
}

// psi_R(x, y) = k(A2, y) * du/dx(A2, y)
double psi_R(double y) {
	return (-12) * (pow((y + 3), 2) + 1) / pow((pow(y, 2) + 10), 2);
}

// psi_L(x, y) = -k(A1, y) * du/dx(A1, y)
double psi_L(double y) {
	return (-8) * (pow((y - 2), 2) + 1) / pow((pow(y, 2) + 5), 2);
}

// psi_T(x, y) = k(x, B2) * du/dy(x, B2)
double psi_T(double x) {
	return (-16) * (pow((x + 4), 2) + 1) / pow((pow(x, 2) + 17), 2);
}

// psi_B(x, y) = -k(x, B2) * du/dy(x, B2)
double psi_B(double x) {
	return (-4) * (pow((x - 1), 2) + 1) / pow((pow(x, 2) + 2), 2);
}

// right-part function of Poisson equation
double F(double x, double y) {
	return 2 * (pow(x,4) + pow(y,4) + 2 * (pow(x,2) + 3) * pow(y,2) + 6 * pow(x,2) + 16*x*y + 5)
	/ pow((1 + pow(x, 2) + pow(y, 2)), 3);
}


double Get(int i, int j){
	return w[j*(M + 1) + i];
}


void CalculateExactSolution(vector<double> &exact) {
	for (int i = 0; i < M + 1; ++i)
		for (int j = 0; j < N + 1; ++j)
			exact[j * (M + 1) + i] = u(A1 + i * h1, B1 + j * h2);

	//cout << "Exact solution of Poisson equation:" << endl;
	//PrintVector(exact);
	//cout << endl;
}


int IsNodeInternalCornerOrSide(int current_node_global_offset_x, int current_node_global_offset_y){

    //corners
    //left bottom corner
    if (current_node_global_offset_x == 0 && current_node_global_offset_y == 0){
        return 2;
    }

    //left top corner
    if (current_node_global_offset_x == 0 && current_node_global_offset_y == N){
        return 4;
    }

    //right bottom corner
    if (current_node_global_offset_x == M && current_node_global_offset_y == 0){
        return 6;
    }

    //right top corner
    if (current_node_global_offset_x == M && current_node_global_offset_y == N){
        return 8;
    }

    //sides
    //left side
    if (current_node_global_offset_y >= 1 && current_node_global_offset_y <= N - 1 &&
        current_node_global_offset_x == 0){
        return 1;
    }

    //right side
    if (current_node_global_offset_y >= 1 && current_node_global_offset_y <= N - 1 &&
        current_node_global_offset_x == M){
        return 3;
    }

    //bottom side
    if (current_node_global_offset_x >= 1 && current_node_global_offset_x <= M - 1 &&
        current_node_global_offset_y == 0){
        return 5;
    }

    //top side
    if ((current_node_global_offset_x >= 1 && current_node_global_offset_x <= M - 1 &&
        current_node_global_offset_y == N)){
        return 7;
    }

    //internal
    if ((current_node_global_offset_x >= 1 && current_node_global_offset_x <= M - 1) &&
        (current_node_global_offset_y >= 1 && current_node_global_offset_y <= N - 1)){
        return 0;
    }

    return -1;
}

double ComputeMagicInnerProductA_iw (int current_node_global_offset_x, int current_node_global_offset_y){

        int glob_x = current_node_global_offset_x;
        int glob_y = current_node_global_offset_y;

        double result = 0.0;

        map <string,bool> neighbours = {
            {"left", true},
            {"right", true},
            {"bottom", true},
            {"top", true}
        };

        double left_neighbour = 0.0, right_neighbour = 0.0, bottom_neighbour = 0.0, top_neighbour = 0.0, this_node = 0.0;
        double left_coeff = 1.0, right_coeff = 1.0, bottom_coeff = 1.0, top_coeff = 1.0, this_coeff = 1.0;

        switch (IsNodeInternalCornerOrSide(glob_x, glob_y)){
            case 2:
                //left bottom corner
                neighbours["left"] = false;
                neighbours["bottom"] = false;
                break;
            case 4:
                //left top corner
                neighbours["left"] = false;
                neighbours["top"] = false;
                break;
            case 6:
                //right bottom corner
                neighbours["right"] = false;
                neighbours["bottom"] = false;
                break;
            case 8:
                //right top corner
                neighbours["right"] = false;
                neighbours["top"] = false;
                break;
            case 1:
                //left side
                neighbours["left"] = false;
                break;
            case 3:
                //right side
                neighbours["right"] = false;
                break;
            case 5:
                //bottom side
                neighbours["bottom"] = false;
                break;
            case 7:
                //top side
                neighbours["top"] = false;
                break;
            case 0:
                //internal
                break;
            default:
                cout << "[ERROR]: Bad global coords compute matrix. Global:" << glob_x << " " << glob_y<<endl;
        }

        if (!neighbours["left"]){
            right_coeff = 2.0;
            left_coeff = 0.0;
        }

        if (!neighbours["right"]){
            left_coeff = 2.0;
            right_coeff = 0.0;
        }

        if (!neighbours["bottom"]){
            top_coeff = 2.0;
            bottom_coeff = 0.0;
        }

        if (!neighbours["top"]){
            bottom_coeff = 2.0;
            top_coeff = 0.0;
        }



        if (neighbours["left"]){
            left_coeff *= -k(x_left + (glob_x - 0.5) * h_x, y_bottom + glob_y * h_y) / pow(h_x, 2);
            left_neighbour = Get(glob_x - 1, glob_y);
        }

        if (neighbours["right"]){
            right_coeff *= -k(x_left + (glob_x + 0.5) * h_x, y_bottom + glob_y * h_y) / pow(h_x, 2);
            right_neighbour = Get(glob_x + 1, glob_y);
        }

        if (neighbours["bottom"]){
            bottom_coeff *= -k(x_left + glob_x * h_x, y_bottom + (glob_y - 0.5) * h_y) / pow(h_y, 2);
            bottom_neighbour = Get(glob_x, glob_y - 1);
        }

        if (neighbours["top"]){
            top_coeff *= -k(x_left + glob_x * h_x, y_bottom + (glob_y + 0.5) * h_y) / pow(h_y, 2);
            top_neighbour = Get(glob_x, glob_y + 1);
        }

        this_coeff = q(x_left + glob_x * h_x, y_bottom + glob_y * h_y) - left_coeff - right_coeff - bottom_coeff - top_coeff;
        this_node = Get(glob_x, glob_y);

        result = left_coeff * left_neighbour +
                 right_coeff * right_neighbour +
                 bottom_coeff * bottom_neighbour +
                 top_coeff * top_neighbour +
                 this_coeff * this_node;

        return result;
    }



    double GetNodeFromB(int current_node_global_offset_x, int current_node_global_offset_y) {

        int glob_x = current_node_global_offset_x;
        int glob_y = current_node_global_offset_y;

        double result = 0.0;

        switch (IsNodeInternalCornerOrSide(glob_x, glob_y)){
            case 2:
                //left bottom corner
                result = F(x_left, y_bottom) + 2.0 / h_x * psi_L(y_bottom) + 2.0 / h_y * psi_B(x_left);
                break;
            case 4:
                //left top corner
                result = F(x_left, y_top) + 2.0 / h_x * psi_L(y_top) + 2.0 / h_y * psi_T(x_left);
                break;
            case 6:
                //right bottom corner
                result = F(x_right, y_bottom) + 2.0 / h_x * psi_R(y_bottom) + 2.0 / h_y * psi_B(x_right);
                break;
            case 8:
                //right top corner
                result = F(x_right, y_top) + 2.0 / h_x * psi_R(y_top) + 2.0 / h_y * psi_T(x_right);
                break;
            case 1:
                //left side
                result = F(x_left, y_bottom + glob_y * h_y) + 2.0 / h_x * psi_L(y_bottom + glob_y * h_y);
                break;
            case 3:
                //right side
                result = F(x_right, y_bottom + glob_y * h_y) + 2.0 / h_x * psi_R(y_bottom + glob_y * h_y);
                break;
            case 5:
                //bottom side
                result = F(x_left + glob_x * h_x, y_bottom) + 2.0 / h_y * psi_B(x_left + glob_x * h_x);
                break;
            case 7:
                //top side
                result = F(x_left + glob_x * h_x, y_top) + 2.0 / h_y * psi_T(x_left + glob_x * h_x);
                break;
            case 0:
                //internal
                result = F(x_left + glob_x * h_x, y_bottom + glob_y * h_y);
                break;
            default:
                cout << "[ERROR]: Bad global coords compute matrix. Global:" << glob_x << " " << glob_y <<endl;

        }

        return result;

    }



double func_inner_product(vector<double> &v1, vector<double> &v2) {
	double prod = 0;

	// rho = 1 in case of interior points
	for (int i = 1; i < M; ++i)
		for (int j = 1; j < N; ++j)
			prod += v1[j * (M + 1) + i] * v2[j * (M + 1) + i];

	// rho = 0.5 in case of left and right sides
	for (int j = 1; j < N; ++j)
		prod += 0.5 * (v1[j * (M + 1)] * v2[j * (M + 1)] + v1[j * (M + 1) + M] * v2[j * (M + 1) + M]);

	// rho = 0.5 in case of top and bottom sides
	for (int i = 1; i < M; ++i)
		prod += 0.5 * (v1[N * (M + 1) + i] * v2[N * (M + 1) + i] + v1[i] * v2[i]);

	// rho = 0.25 in case of corner points
	prod += 0.25 * (v1[0] * v2[0] + v1[M] * v2[M] + v1[N * (M + 1)] * v2[N * (M + 1)]
		+ v1[N * (M + 1) + M] * v2[N * (M + 1) + M]);

	prod *= h1 * h2;
	return prod;
}

double norm(vector<double> &v) {
	return sqrt(func_inner_product(v, v));
}

// solving system by method of least residuals
void LeastResiduals(vector< vector<double> > &A, vector<double> &B) {
	int sz = B.size(), iterations = 0;
	double tau = 0.0, eps = 1e-06;
	vector<double> r(sz, 0.0), Ar(sz, 0.0);

	do {
		// filling residual
		//for (int i = 0; i < sz; ++i)

        for(int i = 0; i < M + 1; ++i){
            for(int j = 0; j < N + 1; ++j){
			//r[i] = inner_product(A[i], w) - B[i];
				r[j * (M + 1) + i] = ComputeMagicInnerProductA_iw(i, j) - GetNodeFromB(i, j);

		old_w = w;

		w = r;

		// filling Ar which is product of matrix A and vector r
		//for (int i = 0; i < sz; ++i)
			//Ar[i] = inner_product(A[i], r);


        for(int i = 0; i < M + 1; ++i){
            for(int j = 0; j < N + 1; ++j){
            	Ar[j * (M + 1) + i] = ComputeMagicInnerProductA_iw(i, j);

		// computing iteration parameter
		tau = func_inner_product(Ar, r) / func_inner_product(Ar, Ar);

		// computing new iteration of w
		for (int i = 0; i < sz; ++i)
			w[i] -= tau * r[i];

		++iterations;

		cout << fabs(tau) * norm(r) - eps << endl;
	} while (fabs(tau) * norm(r) >= eps);

	cout << "Number of iterations: " << iterations << endl;
}

// computing norm of difference between approximate and exact solutions
double Compare(vector<double> &exact) {
	int sz = exact.size();

	// filling difference vector
	vector<double> diff(sz, 0.0);
	for (int i = 0; i < sz; ++i)
		diff[i] = w[i] - exact[i];

	return norm(diff);
}

int main() {
	//cout << "Starting system construction:" << endl;
	vector< vector<double> > A(SIZE, vector<double>(SIZE, 0.0));
	vector<double> B(SIZE, 0.0);

	// construction of the system Aw = B
	ComputeRightPart(B);
	//cout << "[+] Right part constructed" << endl;

	ComputeMatrix(A);
	//cout << "[+] Matrix constructed" << endl;

	// calculation of the exact solution
	vector<double> exact(SIZE, 0.0);
	CalculateExactSolution(exact);
	//cout << "[+] Exact solution calculated" << endl;

	// solving system Aw = B by method of least residuals
	w.resize(SIZE);
	LeastResiduals(A, B);
	cout << "[+] Approximate solution calculated" << endl;

	cout << "Norm of difference between exact and approximate solutions: " << Compare(w, exact) << endl;

	cout << "Done!" << endl;
	return 0;
}