---

/**
 * ============================================================
 *  МЕТОД СЕТОК (Конечно-разностный метод)
 *  Решение уравнения Пуассона на прямоугольной области
 *
 *  Уравнение:  d²u/dx² + d²u/dy² = f(x, y)
 *  Область:    (x, y) ∈ [0, a] × [0, b]
 *  Граничные условия Дирихле: u|_Γ = g(x, y)
 *
 *  Метод решения: Зейделя с верхней релаксацией (SOR)
 * ============================================================
 */

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <fstream>
#include <string>
#include <chrono>

using namespace std;

// ============================================================
//  ОПРЕДЕЛЕНИЕ π (кроссплатформенное решение)
// ============================================================
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

// ============================================================
//  НАСТРОЙКИ ЗАДАЧИ
// ============================================================

// Параметры области
const double A = 1.0;          // Размер по x: [0, A]
const double B = 1.0;          // Размер по y: [0, B]

// Параметры сетки
const int NX = 30;             // Количество узлов по x (включая границы)
const int NY = 30;             // Количество узлов по y (включая границы)

// Параметры итерационного процесса
const double EPS = 1e-6;       // Точность сходимости
const int MAX_ITER = 100000;   // Максимальное число итераций
const double OMEGA = 1.5;      // Параметр релаксации SOR (1 < w < 2)

// ============================================================
//  ФУНКЦИИ ЗАДАЧИ
// ============================================================

/**
 * Правая часть уравнения Пуассона: f(x, y)
 *
 * Тестовый пример:
 *   Точное решение: u(x,y) = sin(pi*x) * sin(pi*y)
 *   Тогда: f(x,y) = -2*pi^2 * sin(pi*x) * sin(pi*y)
 */
double f(double x, double y) {
    return -2.0 * M_PI * M_PI * sin(M_PI * x) * sin(M_PI * y);
}

/**
 * Граничное условие: g(x, y) на границе Γ
 * Для тестового примера: u = 0 на всей границе
 */
double g(double x, double y) {
    return 0.0;
}

/**
 * Точное решение (для проверки точности)
 */
double exact_solution(double x, double y) {
    return sin(M_PI * x) * sin(M_PI * y);
}

// ============================================================
//  КЛАСС: Сеточный решатель
// ============================================================

class GridSolver {
private:
    int nx, ny;           // Число узлов сетки
    double a, b;          // Размеры области
    double hx, hy;        // Шаги сетки
    double omega;         // Параметр релаксации
    double eps;           // Точность
    int max_iter;         // Макс. итераций

    vector<vector<double>> u;       // Решение
    vector<vector<double>> f_val;   // Значения правой части

public:
    // ----------------------------------------------------------
    //  Конструктор
    // ----------------------------------------------------------
    GridSolver(int nx_, int ny_, double a_, double b_,
               double omega_, double eps_, int max_iter_)
        : nx(nx_), ny(ny_), a(a_), b(b_),
          omega(omega_), eps(eps_), max_iter(max_iter_)
    {
        hx = a / (nx - 1);
        hy = b / (ny - 1);

        // Инициализация массивов
        u.assign(nx, vector<double>(ny, 0.0));
        f_val.assign(nx, vector<double>(ny, 0.0));

        // Вычисление значений правой части во внутренних узлах
        for (int i = 1; i < nx - 1; i++) {
            for (int j = 1; j < ny - 1; j++) {
                double x = i * hx;
                double y = j * hy;
                f_val[i][j] = f(x, y);
            }
        }

        // Установка граничных условий
        set_boundary_conditions();
    }

    // ----------------------------------------------------------
    //  Установка граничных условий Дирихле
    // ----------------------------------------------------------
    void set_boundary_conditions() {
        for (int i = 0; i < nx; i++) {
            for (int j = 0; j < ny; j++) {
                if (i == 0 || i == nx - 1 || j == 0 || j == ny - 1) {
                    double x = i * hx;
                    double y = j * hy;
                    u[i][j] = g(x, y);
                }
            }
        }
    }

    // ----------------------------------------------------------
    //  Метод Зейделя с верхней релаксацией (SOR)
    //
    //  Конечно-разностная схема (5-точечный шаблон):
    //
    //  (u[i+1][j] - 2u[i][j] + u[i-1][j]) / hx^2 +
    //  (u[i][j+1] - 2u[i][j] + u[i][j-1]) / hy^2 = f[i][j]
    // ----------------------------------------------------------
    int solve() {
        double cx = 1.0 / (hx * hx);
        double cy = 1.0 / (hy * hy);
        double denominator = 2.0 * cx + 2.0 * cy;

        int iter = 0;
        double max_diff;

        cout << "Запуск итерационного процесса (SOR, w = "
             << omega << ")..." << endl;

        for (iter = 1; iter <= max_iter; iter++) {
            max_diff = 0.0;

            for (int i = 1; i < nx - 1; i++) {
                for (int j = 1; j < ny - 1; j++) {
                    // Новое значение по формуле Зейделя
                    double u_new = (
                        (u[i + 1][j] + u[i - 1][j]) * cx +
                        (u[i][j + 1] + u[i][j - 1]) * cy -
                        f_val[i][j]
                    ) / denominator;

                    // SOR: релаксация
                    u_new = omega * u_new + (1.0 - omega) * u[i][j];

                    // Отслеживание максимальной невязки
                    double diff = fabs(u_new - u[i][j]);
                    if (diff > max_diff) {
                        max_diff = diff;
                    }

                    u[i][j] = u_new;
                }
            }

            // Вывод прогресса каждые 500 итераций
            if (iter % 500 == 0) {
                cout << "  Итерация " << setw(6) << iter
                     << " | max|du| = " << scientific
                     << setprecision(4) << max_diff << endl;
            }

            // Проверка сходимости
            if (max_diff < eps) {
                break;
            }
        }

        return iter;
    }

    // ----------------------------------------------------------
    //  Вычисление погрешности
    // ----------------------------------------------------------
    double compute_error() {
        double max_error = 0.0;

        for (int i = 0; i < nx; i++) {
            for (int j = 0; j < ny; j++) {
                double x = i * hx;
                double y = j * hy;
                double error = fabs(u[i][j] - exact_solution(x, y));
                if (error > max_error) {
                    max_error = error;
                }
            }
        }

        return max_error;
    }

    // ----------------------------------------------------------
    //  Вывод решения в консоль (тепловая карта)
    // ----------------------------------------------------------
    void print_heatmap() {
        cout << "\n=== ТЕПЛОВАЯ КАРТА РЕШЕНИЯ ===" << endl;
        cout << "(символы: ' ' < '.' < ':' < '-' < '=' < '+' < '*' < '#')"
             << endl << endl;

        const string chars = " .:-=+*#%@";
        int levels = (int)chars.size();

        // Находим min/max для масштабирования
        double u_min = 1e18, u_max = -1e18;
        for (int i = 0; i < nx; i++) {
            for (int j = 0; j < ny; j++) {
                if (u[i][j] < u_min) u_min = u[i][j];
                if (u[i][j] > u_max) u_max = u[i][j];
            }
        }

        double range = u_max - u_min;
        if (range < 1e-15) range = 1.0;

        // Вывод (y сверху вниз для корректного отображения)
        for (int j = ny - 1; j >= 0; j--) {
            for (int i = 0; i < nx; i++) {
                int level = (int)((u[i][j] - u_min) / range * (levels - 1));
                if (level < 0) level = 0;
                if (level >= levels) level = levels - 1;
                cout << chars[level];
            }
            if (j == ny - 1) cout << "  y = " << fixed << setprecision(2) << b;
            if (j == 0)      cout << "  y = 0";
            cout << endl;
        }

        cout << "x = 0"
             << string(nx > 10 ? nx - 10 : 0, ' ')
             << "x = " << fixed << setprecision(2) << a << endl;
        cout << "\nМинимум: " << scientific << u_min
             << "  Максимум: " << u_max << endl;
    }

    // ----------------------------------------------------------
    //  Вывод числовой таблицы
    // ----------------------------------------------------------
    void print_table(int step = 1) {
        cout << "\n=== ЧИСЛОВАЯ ТАБЛИЦА РЕШЕНИЯ ===" << endl;
        cout << setw(8) << "x\\y";

        for (int j = 0; j < ny; j += step) {
            cout << setw(10) << fixed << setprecision(3) << j * hy;
        }
        cout << endl;

        int cols = (ny - 1) / step + 1;
        cout << string(8 + 10 * cols, '-') << endl;

        for (int i = 0; i < nx; i += step) {
            cout << setw(8) << fixed << setprecision(3) << i * hx;
            for (int j = 0; j < ny; j += step) {
                cout << setw(10) << fixed << setprecision(5) << u[i][j];
            }
            cout << endl;
        }
    }

    // ----------------------------------------------------------
    //  Экспорт в файл
    // ----------------------------------------------------------
    void export_to_file(const string& filename) {
        ofstream file(filename);
        if (!file.is_open()) {
            cerr << "Ошибка: не удалось открыть файл " << filename << endl;
            return;
        }

        file << "# x  y  u(x,y)  u_exact  error" << endl;
        for (int i = 0; i < nx; i++) {
            for (int j = 0; j < ny; j++) {
                double x = i * hx;
                double y = j * hy;
                file << fixed << setprecision(8)
                     << x << "  " << y << "  "
                     << u[i][j] << "  "
                     << exact_solution(x, y) << "  "
                     << fabs(u[i][j] - exact_solution(x, y))
                     << endl;
            }
            file << endl;
        }

        file.close();
        cout << "Результаты экспортированы в файл: " << filename << endl;
    }

    // ----------------------------------------------------------
    //  Информация о сетке
    // ----------------------------------------------------------
    void print_info() {
        cout << "\n+--------------------------------------------------+" << endl;
        cout << "|        МЕТОД СЕТОК - УРАВНЕНИЕ ПУАССОНА          |" << endl;
        cout << "+--------------------------------------------------+" << endl;
        cout << "|  Область:    [0, " << A << "] x [0, " << B << "]"
             << string(20, ' ') << "|" << endl;
        cout << "|  Сетка:      " << nx << " x " << ny << " узлов"
             << string(24, ' ') << "|" << endl;
        cout << "|  Шаг hx:     " << fixed << setprecision(6) << hx
             << string(18, ' ') << "|" << endl;
        cout << "|  Шаг hy:     " << fixed << setprecision(6) << hy
             << string(18, ' ') << "|" << endl;
        cout << "|  Точность:   " << scientific << eps
             << string(18, ' ') << "|" << endl;
        cout << "|  Параметр w: " << fixed << setprecision(2) << omega
             << string(23, ' ') << "|" << endl;
        cout << "+--------------------------------------------------+" << endl;
    }
};

// ============================================================
//  ИССЛЕДОВАНИЕ СХОДИМОСТИ ПО РАЗМЕРУ СЕТКИ
// ============================================================

void convergence_study() {
    cout << "\n+--------------------------------------------------+" << endl;
    cout << "|       ИССЛЕДОВАНИЕ СХОДИМОСТИ ПО СЕТКЕ           |" << endl;
    cout << "+--------------------------------------------------+" << endl;
    cout << "|  N    |  hx       |  Итераций |  Max погрешность |" << endl;
    cout << "+-------+-----------+-----------+------------------+" << endl;

    vector<int> grid_sizes = {10, 15, 20, 30, 40, 50};

    for (int n : grid_sizes) {
        GridSolver solver(n, n, A, B, OMEGA, EPS, MAX_ITER);
        int iters = solver.solve();
        double error = solver.compute_error();

        cout << "|  " << setw(3) << n << "  | "
             << fixed << setprecision(5) << A / (n - 1) << "  | "
             << setw(7) << iters << "   | "
             << scientific << setprecision(4) << error
             << "    |" << endl;
    }

    cout << "+--------------------------------------------------+" << endl;
}

// ============================================================
//  ГЛАВНАЯ ФУНКЦИЯ
// ============================================================

int main() {
    setlocale(LC_ALL, "Russian");

    cout << "=====================================================" << endl;
    cout << "   ЧИСЛЕННЫЕ МЕТОДЫ: МЕТОД СЕТОК (ВАРИАНТ 2)" << endl;
    cout << "   Решение уравнения Пуассона методом SOR" << endl;
    cout << "=====================================================" << endl;

    // ---- Основной расчёт ----
    GridSolver solver(NX, NY, A, B, OMEGA, EPS, MAX_ITER);
    solver.print_info();

    // Замер времени
    auto start = chrono::high_resolution_clock::now();

    int iterations = solver.solve();

    auto end = chrono::high_resolution_clock::now();
    double elapsed = chrono::duration<double, milli>(end - start).count();

    // ---- Результаты ----
    double max_error = solver.compute_error();

    cout << "\n+--------------------------------------------------+" << endl;
    cout << "|  РЕЗУЛЬТАТЫ                                      |" << endl;
    cout << "+--------------------------------------------------+" << endl;
    cout << "|  Итераций выполнено:  " << setw(8) << iterations
         << string(18, ' ') << "|" << endl;
    cout << "|  Время вычислений:    " << fixed << setprecision(2)
         << elapsed << " мс" << string(16, ' ') << "|" << endl;
    cout << "|  Макс. погрешность:   " << scientific << setprecision(6)
         << max_error << string(8, ' ') << "|" << endl;
    cout << "|  Порядок точности:    O(h^2) = O("
         << fixed << setprecision(6) << (A / (NX - 1)) * (A / (NX - 1))
         << ")" << string(6, ' ') << "|" << endl;
    cout << "+--------------------------------------------------+" << endl;

    // ---- Визуализация ----
    solver.print_heatmap();

    // Числовая таблица (с шагом 5 для читаемости)
    solver.print_table(5);

    // ---- Экспорт ----
    solver.export_to_file("grid_solution.dat");

    // ---- Исследование сходимости ----
    convergence_study();

    cout << "\nПрограмма завершена." << endl;

    return 0;
}