---

/**
 * ============================================================
 *  МЕТОД СЕТОК (Конечно-разностный метод)
 *  Решение уравнения Пуассона на прямоугольной области
 *
 *  Уравнение:  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 =