УМФ

#include <iostream>
#include <fstream>
#include <functional>
#include <vector>
#include <string>
#include <tuple>

using std::cout;
using std::ifstream;
using std::ofstream;
using std::endl;
using std::function;
using std::vector;
using std::string;
using std::ostream;
using std::istream;
using std::swap;
using std::tuple;
using std::make_tuple;
using std::tie;

#define DIM_LOCAL 4

//--------------------------------------------------------------------------------
// Структура для хранение тестировочных данных
// U        - истиннове решение (x, y, t)
// f        - функция правой части (x, y, t)
// sigma    - сигма из уравнения
// hee      - хи из уравнения
// lambda   - лямбда из уравнения
//--------------------------------------------------------------------------------
struct TestData {
    function<double(double, double, double)> u;
    function<double(double, double, double)> f;
    double sigma;
    double hee;
    double lambda;
};

//--------------------------------------------------------------------------------
// Структура для хранения простейших первых граничных условий,
// прямоугольная область, на пересечениях границ функции должны быть равны
// leftXCondition   - левая граница области(x = x0)
// rightXCondition  - правая граница области(x = x_max)
// leftYCondition   - нижняя граница области(y = y0)
// rightYCondition  - верхняя граница области(y = y_max)
//--------------------------------------------------------------------------------
struct Temp1stBoundCondition {
    function<double(double, double, double)> leftXCondition;
    function<double(double, double, double)> rightXCondition;
    function<double(double, double, double)> leftYCondition;
    function<double(double, double, double)> rightYCondition;
};

//--------------------------------------------------------------------------------
// Функция загрузки и конструирования сетки
// Формат данных в потоке: начальная точка, конечная точка,
//  количество точек с начальной и конечной, коэффициент растяжения/сжатия
// in       - входной поток
// grid     - выход, координаты узлов сетки
//--------------------------------------------------------------------------------
void uploadGrid(istream & in, vector<double> & grid) {
    double start, end, n, k;
    double h;
    in >> start >> end >> n >> k;
    n = n - 1;
    if (!grid.empty()) grid.clear();

    grid.push_back(start);
    if (k != 1.) {
        h = (end - start) * (1 - k) / (1 - pow(k, n));
        for (int i = 1; i < n; i++) {
            grid.push_back(grid.back() + h);
            h *= k;
        }
    }
    else {
        h = (end - start) / (n);
        for (int i = 1; i < n; i++) {
            grid.push_back(grid.back() + h);
        }
    }
    grid.push_back(end);
}

//--------------------------------------------------------------------------------
// Постройка сеток для двумерной нестационарной задачи
// xGrid        - выход, сетка по иксу
// yGrid        - выход, сетка по игреку
// tGrid        - выход, сетка по времени
// xGridName    - имя файла с параметрами сетки по иксу
// yGridName    - имя файла с параметрами сетки по игреку
// tGridName    - имя файла с параметрами сетки по времени
//--------------------------------------------------------------------------------
void constructGrids(vector<double> & xGrid, vector<double> & yGrid, vector<double> & tGrid,
                    const string & xGridName = "grid_x.txt", const string & yGridName = "grid_y.txt", const string & tGridName = "grid_t.txt") {
    ifstream in(xGridName);
    uploadGrid(in, xGrid);
    in.close();
    in.open(yGridName);
    uploadGrid(in, yGrid);
    in.close();
    in.open(tGridName);
    uploadGrid(in, tGrid);
    in.close();
}

//--------------------------------------------------------------------------------
// Ядро локальной матрицы Жёсткости(с множителем l/6*hy/hx)
//--------------------------------------------------------------------------------
double kernGLocalYX[DIM_LOCAL][DIM_LOCAL] = {   { 2, -2, 1, -1},
                                                { -2, 2, -1, 1},
                                                { 1, -1, 2, -2},
                                                { -1, 1, -2, 2} };

//--------------------------------------------------------------------------------
// Ядро локальной матрицы Жёсткости(с множителем l/6*hx/hy)
//--------------------------------------------------------------------------------
double kernGLocalXY[DIM_LOCAL][DIM_LOCAL] = {   { 2, 1, -2, -1 },
                                                { 1, 2, -1, -2 },
                                                { -2, -1, 2, 1 },
                                                { -1, -2, 1, 2 } };

//--------------------------------------------------------------------------------
// Ядро локальной матрицы Масс(Ядро локальной матрицы C)
//--------------------------------------------------------------------------------
double kernMLocal[DIM_LOCAL][DIM_LOCAL] = { { 4, 2, 2, 1 },
                                            { 2, 4, 1, 2 },
                                            { 2, 1, 4, 2 },
                                            { 1, 2, 2, 4 } };

//--------------------------------------------------------------------------------
// Получение компоненты локального вектора b
// index    - элемент вектора b(счёт с нуля)
// f        - функция правой части (x, y, t)
// x1       - начальная координата по иксу конечного элемента
// x2       - конечная координата по иксу конечного элемента
// y1       - начальная координата по игреку конечного элемента
// y2       - конечная координата по игреку конечного элемента
// t        - значение времени текущего временного слоя
//--------------------------------------------------------------------------------
double getBLocal(const size_t index, const function<double(double, double, double)> & f,
                const double x1, const double x2, const double y1, const double y2, const double t) {
    double result = 0;

    result += kernMLocal[index][0] * f(x1, y1, t);
    result += kernMLocal[index][1] * f(x2, y1, t);
    result += kernMLocal[index][2] * f(x1, y2, t);
    result += kernMLocal[index][3] * f(x2, y2, t);

    return result;// *(x2 - x1) * (y2 - y1) / 36 * (x2 - x1) * (y2 - y1) / 36;
}

//--------------------------------------------------------------------------------
// Вставляет локальную матрицу в глобальную
// local    - локальная матрица
// global   - глобальная матрица
// i1       - первый глобальный индекс вставки
// i2       - второй глобальный индекс вставки
// i3       - третий глобальный индекс вставки
// i4       - четвёртый глобальный индекс вставки
//--------------------------------------------------------------------------------
void insertLocalToGlobal(const double (*local)[DIM_LOCAL], double ** global, const double mult,
    const size_t i1, const size_t i2, const size_t i3, const size_t i4) {

    vector<size_t> indexes = { i1, i2, i3, i4 };
    for (int i = 0; i < DIM_LOCAL; i++) {
        for (int j = 0; j < DIM_LOCAL; j++) {
            global[indexes[i]][indexes[j]] += mult * local[i][j];
        }
    }
}

//--------------------------------------------------------------------------------
// Вставляет локальную матрицу в глобальную
// local    - локальная матрица
// global   - глобальная матрица
// i1       - первый глобальный индекс вставки
// i2       - второй глобальный индекс вставки
// i3       - третий глобальный индекс вставки
// i4       - четвёртый глобальный индекс вставки
//--------------------------------------------------------------------------------
void insertLocalToGlobal(const double * local, double * global, const double mult,
    const size_t i1, const size_t i2, const size_t i3, const size_t i4) {

    global[i1] += mult * local[0];
    global[i2] += mult * local[1];
    global[i3] += mult * local[2];
    global[i4] += mult * local[3];
}

//--------------------------------------------------------------------------------
// Вспомогательная функция для вычисление глобального индекса первого узла по
// номеру конечного элемента
// numberFE - номер конечного элемента(счёт с 0)
// sizeX    - количество узлов по иксу
//--------------------------------------------------------------------------------
inline size_t globalStartIndex(const size_t numberFE, const size_t sizeX) {
    return numberFE + static_cast<size_t>(numberFE / (sizeX - 1));
}

//--------------------------------------------------------------------------------
// Вспомогательная функция для вычисление индекса икса по глобальному индексу узла
// numberFE - номер конечного элемента(счёт с 0)
// sizeX    - количество узлов по иксу
//--------------------------------------------------------------------------------
inline size_t indexXbyGlobalIndex(const size_t numberFE, const size_t sizeX) {
    return globalStartIndex(numberFE,sizeX) % sizeX;
}

//--------------------------------------------------------------------------------
// Вспомогательная функция для вычисление индекса игрека по глобальному индексу
// узла
// numberFE - номер конечного элемента(счёт с 0)
// sizeX    - количество узлов по иксу
//--------------------------------------------------------------------------------
inline size_t indexYbyGlobalIndex(const size_t numberFE, const size_t sizeX) {
    return static_cast<size_t>(globalStartIndex(numberFE, sizeX) / sizeX);
}

//--------------------------------------------------------------------------------
// Сборка глобальной матрицы Масс(гамма)
// globalM  - глобальная матрица Масс
// xGrid    - сетка по иксу
// yGrid    - сетка по игреку
// gamma    - гамма(множитель)
//--------------------------------------------------------------------------------
void constructGlobalM(double ** globalM, const vector<double> & xGrid,
    const vector<double> & yGrid, const double gamma) {

    size_t countFE = (xGrid.size() - 1) * (yGrid.size() - 1);
    double dx, dy;
    double xStart, xEnd;
    double yStart, yEnd;
    double globalStart;

    for (int elemNum = 0; elemNum < countFE; elemNum++) {
        globalStart = globalStartIndex(elemNum, xGrid.size());
        xStart = xGrid[indexXbyGlobalIndex(elemNum, xGrid.size())];
        xEnd = xGrid[indexXbyGlobalIndex(elemNum, xGrid.size()) + 1];
        yStart = yGrid[indexYbyGlobalIndex(elemNum, xGrid.size())];
        yEnd = yGrid[indexYbyGlobalIndex(elemNum, xGrid.size()) + 1];

        dx = xEnd - xStart;
        dy = yEnd - yStart;

        insertLocalToGlobal(kernMLocal, globalM, gamma * dx * dy / 36,
            globalStart, globalStart + 1, globalStart + xGrid.size(), globalStart + xGrid.size() + 1);
    }
}

//--------------------------------------------------------------------------------
// Сборка глобальный матрицы Жёсткости(лямбда)
// globalM  - глобальная матрица Масс
// xGrid    - сетка по иксу
// yGrid    - сетка по игреку
// lambda   - лямбда(множитель)
//--------------------------------------------------------------------------------
void constructGlobalG(double ** globalG, const vector<double> & xGrid,
    const vector<double> & yGrid, const double lambda) {

    size_t countFE = (xGrid.size() - 1) * (yGrid.size() - 1);
    double dx, dy;
    double xStart, xEnd;
    double yStart, yEnd;
    double globalStart;

    for (int elemNum = 0; elemNum < countFE; elemNum++) {
        globalStart = globalStartIndex(elemNum, xGrid.size());
        xStart = xGrid[indexXbyGlobalIndex(elemNum, xGrid.size())];
        xEnd = xGrid[indexXbyGlobalIndex(elemNum, xGrid.size()) + 1];
        yStart = yGrid[indexYbyGlobalIndex(elemNum, xGrid.size())];
        yEnd = yGrid[indexYbyGlobalIndex(elemNum, xGrid.size()) + 1];

        dx = xEnd - xStart;
        dy = yEnd - yStart;

        insertLocalToGlobal(kernGLocalYX, globalG, lambda * dy / 6. / dx,
            globalStart, globalStart + 1, globalStart + xGrid.size(), globalStart + xGrid.size() + 1);
        insertLocalToGlobal(kernGLocalXY, globalG, lambda * dx / 6. / dy,
            globalStart, globalStart + 1, globalStart + xGrid.size(), globalStart + xGrid.size() + 1);
    }
}

//--------------------------------------------------------------------------------
// Сборка глобального вектора правой части
// globalB  - глобальный вектор B
// xGrid    - сетка по иксу
// yGrid    - сетка по игреку
// f        - функция правой части
// t        - текущий временной слой
//--------------------------------------------------------------------------------
void constructGlobalB(double * globalB, const vector<double> & xGrid,
    const vector<double> & yGrid, const function<double(double,double,double)> f,
    const double t) {

    size_t countFE = (xGrid.size() - 1) * (yGrid.size() - 1);
    double dx, dy;
    double xStart, xEnd;
    double yStart, yEnd;
    double globalStart;

    double localRightPart[DIM_LOCAL];

    for (int elemNum = 0; elemNum < countFE; elemNum++) {
        globalStart = globalStartIndex(elemNum, xGrid.size());
        xStart = xGrid[indexXbyGlobalIndex(elemNum, xGrid.size())];
        xEnd = xGrid[indexXbyGlobalIndex(elemNum, xGrid.size()) + 1];
        yStart = yGrid[indexYbyGlobalIndex(elemNum, xGrid.size())];
        yEnd = yGrid[indexYbyGlobalIndex(elemNum, xGrid.size()) + 1];

        dx = xEnd - xStart;
        dy = yEnd - yStart;

        for (int i = 0; i < DIM_LOCAL; i++)
            localRightPart[i] = getBLocal(i, f, xStart, xEnd, yStart, yEnd, t);

        insertLocalToGlobal(localRightPart, globalB, dx * dy / 36.,
            globalStart, globalStart + 1, globalStart + xGrid.size(), globalStart + xGrid.size() + 1);
    }
}

//--------------------------------------------------------------------------------
// Решатель СЛАУ методом Гаусс.
// globalA  - глобальная матрица A
// globalB  - вектор правой части(вход), результат решения(выход)
// n        - размернойсть квадратной матрицы и вектора правой части
//--------------------------------------------------------------------------------
void solveSLAE(const double * const * globalA, double *globalB, size_t n) {
    double tmp;
    double **A = new double *[n];
    for (int i = 0; i < n; i++) {
        A[i] = new double[n];
        for (int j = 0; j < n; j++)
            A[i][j] = globalA[i][j];
    }

    for (int i = 0; i < n; i++) {
        tmp = A[i][i];
        for (int j = i; j < n; j++)
            A[i][j] /= tmp;
        globalB[i] /= tmp;
        for (int j = i + 1; j < n; j++) {
            tmp = A[j][i];
            for (int k = i; k < n; k++)
                A[j][k] -= tmp * A[i][k];
            globalB[j] -= tmp * globalB[i];
        }
        for (int j = i - 1; j >= 0; j--) {
            tmp = A[j][i];
            for (int k = i; k < n; k++)
                A[j][k] -= tmp * A[i][k];
            globalB[j] -= tmp * globalB[i];
        }
    }

    for (int i = 0; i < n; i++)
        delete[] A[i];
    delete[] A;
}

//--------------------------------------------------------------------------------
// Применяем первый краевые условия
// globalA      - глобальная матрица для СЛАУ
// globalB      - глобальный вектор правой части для СЛАУ
// xGrid        - сетка по иксу
// yGrid        - сетка по игреку
// condition    - структура с краевыми условиями первого рода
// t            - текущий временной слой
//--------------------------------------------------------------------------------
void apply1stCondition(double **globalA, double *globalB,
    const vector<double> & xGrid, const vector<double> & yGrid,
    const Temp1stBoundCondition & condition, const double t) {

    size_t countNode = (xGrid.size()) * (yGrid.size());

    // Нижняя граница. y = y_min, x = [x_min, x_max]
    for (int i = 0; i < xGrid.size(); i++) {
        for (int j = 0; j < countNode; j++) {
            globalA[i][j] = 0;
        }
        globalA[i][i] = 1.;

        globalB[i] = condition.leftYCondition(xGrid[i], yGrid[0], t);
    }
    // Верхняя граница. y = y_max, x = [x_min, x_max]
    for (int xIndex = 0, curElem = countNode - xGrid.size(); curElem < countNode; xIndex++, curElem++) {
        for (int j = 0; j < countNode; j++) {
            globalA[curElem][j] = 0;
        }
        globalA[curElem][curElem] = 1.;

        globalB[curElem] = condition.rightYCondition(xGrid[xIndex], yGrid.back(), t);
    }
    // Левая граница. x = x_min, y = [y_min, y_max]
    for (int yIndex = 0, curElem = 0; curElem < countNode; yIndex++, curElem += xGrid.size()) {
        for (int j = 0; j < countNode; j++) {
            globalA[curElem][j] = 0;
        }
        globalA[curElem][curElem] = 1.;

        globalB[curElem] = condition.leftXCondition(xGrid[0], yGrid[yIndex], t);
    }
    // Правая граница. x = x_max, y = [y_min, y_max]
    for (int yIndex = 0, curElem = xGrid.size() - 1; curElem < countNode; yIndex++, curElem += xGrid.size()) {
        for (int j = 0; j < countNode; j++) {
            globalA[curElem][j] = 0;
        }
        globalA[curElem][curElem] = 1.;

        globalB[curElem] = condition.rightXCondition(xGrid.back(), yGrid[yIndex], t);
    }
}

//--------------------------------------------------------------------------------
// Постройка весового вектора по истинному решению
// qi       - глобальный вектор весов
// ti       - значение времени на временном слое
// xGrid    - сетка по иксу
// yGrid    - сетка по игреку
// u        - истинное решение
//--------------------------------------------------------------------------------
void constructQiSyntetic(double * qi, const double ti, const vector<double> & xGrid,
    const vector<double> & yGrid, const function<double(double, double, double)> & u) {

    for (int curY = 0, i = 0; curY < yGrid.size(); curY++) {
        for (int curX = 0; curX < xGrid.size(); curX++, i++) {
            qi[i] = u(xGrid[curX], yGrid[curY], ti);
        }
    }
}

//--------------------------------------------------------------------------------
// Сборка начальных слоёв для трёхточечной схемы
// q_2      - выход, вектор весов соответствующий t = t_0
// q_1      - выход, вектор весов соответствующий t = t_1
// t_2      - значение времени t_0
// t_1      - значение времени t_1
// xGrid    - сетка по иксу
// yGrid    - сетка по игреку
// u        - истинное решение
//--------------------------------------------------------------------------------
void constructStartCondition(double * q_2, double * q_1, const double t_2, const double t_1,
    const vector<double> & xGrid, const vector<double> & yGrid,
    const function<double(double, double, double)> & u) {

    constructQiSyntetic(q_2, t_2, xGrid, yGrid, u);
    constructQiSyntetic(q_1, t_1, xGrid, yGrid, u);

}

//--------------------------------------------------------------------------------
// Правило сборки глобальное матрицы СЛАУ для гиперболического уравнения со схемой
// Кранка-Николсона(трёхточечная неявная по времени)
// globalG      - глобальная марица жёсткости уже умноженная на лямбду
// globalMhi    - глобальная матрица масс умноженная на хи
// globalMsig   - глобальная матрица масс умноженная на сигму
// globalA      - глобальная матрица СЛАУ
// dt           - шаг по времени
// dim          - количество узлов в глобальной сетки(размерность весового вектора)
//--------------------------------------------------------------------------------
void constructGlobalA(const double * const * globalG, const double * const * globalMhi,
    const double * const * globalMsig, double * const * globalA, const double dt, const size_t dim) {

    for (int row = 0; row < dim; row++) {
        for (int col = 0; col < dim; col++) {
            globalA[row][col] = 1. / dt / dt * globalMhi[row][col]
                + 1. / 2. / dt * globalMsig[row][col] + 1. / 2. * globalG[row][col];
        }
    }

}

//--------------------------------------------------------------------------------
// TODO: и это дозаполнять
//--------------------------------------------------------------------------------
void constructGlobalRightPart(const double * globalB, const double * globalB_2,
    const double * const * globalMhi, const double * const * globalMsig, const double *const * globalG,
    const double * qj_1, const double * qj_2, const double dt, double * globalRight,
    const size_t dim) {

    double temp;

    for (int curLine = 0; curLine < dim; curLine++) {
        temp = 0;
        // + 1/2 bj
        globalRight[curLine] = 1. / 2. * globalB[curLine];
        // + 1/2 bj_2
        globalRight[curLine] += 1. / 2. * globalB_2[curLine];
        // + 2. / dt^2 * Mhi * qj-1
        for (int sumIndex = 0; sumIndex < dim; sumIndex++) {
            temp += globalMhi[curLine][sumIndex] * qj_1[sumIndex];
        }
        globalRight[curLine] += 2. / dt / dt * temp;
        temp = 0;
        // - 1./dt^2 * Mhi * qj-2
        for (int sumIndex = 0; sumIndex < dim; sumIndex++) {
            temp += globalMhi[curLine][sumIndex] * qj_2[sumIndex];
        }
        globalRight[curLine] -= 1. / dt / dt * temp;
        temp = 0;
        // + 1/2dt * Msig * qj-2
        for (int sumIndex = 0; sumIndex < dim; sumIndex++) {
            temp += globalMsig[curLine][sumIndex] * qj_2[sumIndex];
        }
        globalRight[curLine] += 1. / 2. / dt * temp;
        temp = 0;
        // - 1/2 * G qj-2
        for (int sumIndex = 0; sumIndex < dim; sumIndex++) {
            temp += globalG[curLine][sumIndex] * qj_2[sumIndex];
        }
        globalRight[curLine] -= 1. / 2. * temp;
    }

}


//--------------------------------------------------------------------------------
// Вывод в файл правой части
// filename     - имя файла
// rightPart    - вектор правой части
// dim          - размерность вектора
//--------------------------------------------------------------------------------
void outputRightPart(const string & fileName, const double * rightPart, const size_t dim) {
    ofstream out(fileName);
    for (int i = 0; i < dim; i++)
        out << rightPart[i] << endl;
    out.close();
}

//--------------------------------------------------------------------------------
// TODO: Заполнить
//--------------------------------------------------------------------------------
void calcErr(const string & fileTruth, const string & curFile, const string & outFile) {
    ifstream inTruth(fileTruth);
    ifstream inCur(curFile);
    ofstream out(outFile);

    while (!inTruth.eof()) {
        out << inTruth.get() - inCur.get() << endl;
    }

    out.close();
    inCur.close();
    inTruth.close();
}

//--------------------------------------------------------------------------------
// TODO: Заполнить
//--------------------------------------------------------------------------------
void swapVectors(double * v1, double * v2, const size_t dim) {
    for (int i = 0; i < dim; i++)
        swap(v1[i], v2[i]);
}

//--------------------------------------------------------------------------------
// TODO: Заполнить и это не забыть
//--------------------------------------------------------------------------------
void clearVector(double * v, const size_t dim) {
    for (int i = 0; i < dim; i++)
        v[i] = 0;
}

//--------------------------------------------------------------------------------
// TODO: Заполнить и это не забыть
//--------------------------------------------------------------------------------
void testing(const TestData & test, const string & testName, const string & xFileName,
             const string & yFileName, const string & tFileName,
             const Temp1stBoundCondition & condition) {

    vector<double> xGrid, yGrid, tGrid;
    constructGrids(xGrid, yGrid, tGrid, xFileName, yFileName, tFileName);

    size_t countFE = (xGrid.size() - 1) * (yGrid.size() - 1);
    size_t countNode = xGrid.size() * yGrid.size();

    double ** globalMhi = new double*[countNode];
    for (int i = 0; i < countNode; i++)
        globalMhi[i] = new double[countNode]();

    double ** globalMsig = new double*[countNode];
    for (int i = 0; i < countNode; i++)
        globalMsig[i] = new double[countNode]();

    double ** globalG = new double*[countNode];
    for (int i = 0; i < countNode; i++)
        globalG[i] = new double[countNode]();

    double ** globalA = new double*[countNode];
    for (int i = 0; i < countNode; i++)
        globalA[i] = new double[countNode]();

    double * curGlobalB = new double[countNode];
    double * curGlobalB_2 = new double[countNode];
    double * curQ = new double[countNode];
    double * q_1 = new double[countNode];
    double * q_2 = new double[countNode];
    double * globalRight = new double[countNode];


    constructStartCondition(q_2, q_1, tGrid[0], tGrid[1], xGrid, yGrid, test.u);
    constructGlobalM(globalMhi, xGrid, yGrid, test.hee);
    constructGlobalM(globalMsig, xGrid, yGrid, test.sigma);
    constructGlobalG(globalG, xGrid, yGrid, test.lambda);

    outputRightPart(testName + "\\t_" + std::to_string(tGrid[0]) + ".txt", q_2, countNode);
    outputRightPart(testName + "\\t_" + std::to_string(tGrid[1]) + ".txt", q_1, countNode);

    double dt;
    for (int curT = 2; curT < tGrid.size(); curT++) {
        dt = tGrid[curT] - tGrid[curT - 1];
        constructGlobalA(globalG, globalMhi, globalMsig, globalA, dt, countNode);
        clearVector(curGlobalB, countNode);
        clearVector(curGlobalB_2, countNode);

        constructGlobalB(curGlobalB, xGrid, yGrid, test.f, tGrid[curT]);
        constructGlobalB(curGlobalB_2, xGrid, yGrid, test.f, tGrid[curT - 2]);
        constructGlobalRightPart(curGlobalB, curGlobalB_2, globalMhi, globalMsig, globalG, q_1, q_2, dt, curQ, countNode);
        apply1stCondition(globalA, curQ, xGrid, yGrid, condition, tGrid[curT]);
        solveSLAE(globalA, curQ, countNode);

        swapVectors(q_2, q_1, countNode);
        swapVectors(q_1, curQ, countNode);
        constructQiSyntetic(curQ, tGrid[curT], xGrid, yGrid, test.u);

        outputRightPart(testName + "\\t_" + std::to_string(tGrid[curT]) + ".txt", q_1, countNode);
        outputRightPart(testName + "\\truth_t_" + std::to_string(tGrid[curT]) + ".txt", curQ, countNode);
        for (int i = 0; i < countNode; i++) {
            curQ[i] -= q_1[i];
        }
        outputRightPart(testName + "\\err_t_" + std::to_string(tGrid[curT]) + ".txt", curQ, countNode);
        outputRightPart(testName + "\\last_error.txt", curQ, countNode);
    }

    double sum = 0.0;
    for (int i = 0; i < countNode; i++) {
        sum += curQ[i];
    }
    sum = sum / countNode;
    std::cout << "Погрешность: " << sum << std::endl;

    delete[] globalRight;
    delete[] q_2;
    delete[] q_1;
    delete[] curQ;
    delete[] curGlobalB_2;
    delete[] curGlobalB;

    for (int i = 0; i < countNode; i++)
        delete[] globalA[i];
    delete[] globalA;

    for (int i = 0; i < countNode; i++)
        delete[] globalG[i];
    delete[] globalG;

    for (int i = 0; i < countNode; i++)
        delete[] globalMsig[i];
    delete[] globalMsig;

    for (int i = 0; i < countNode; i++)
        delete[] globalMhi[i];
    delete[] globalMhi;

}


//--------------------------------------------------------------------------------
//
//--------------------------------------------------------------------------------

int main() {

    // TODO: запустить тесты, отладить.
    TestData test1;
    test1.hee = 1;
    test1.lambda = 6;//1;
    test1.sigma = 1;
    test1.u = [](double x, double y, double t) -> double {
        return x * y*(t - 1);
    };
    test1.f = [](double x, double y, double t) -> double {
        return x * y;
    };
    Temp1stBoundCondition condition1;
    condition1.leftXCondition = [](double x, double y, double t) {
        return  x * y * (t - 1);
    };
    condition1.leftYCondition = [](double x, double y, double t) {
        return  x * y * (t - 1);
    };
    condition1.rightXCondition = [](double x, double y, double t) {
        return x * y * (t - 1);
    };
    condition1.rightYCondition = [](double x, double y, double t) {
        return  x * y * (t - 1);
    };
    //testing(test1, "C:/Users/Alex/Downloads/SuperMFE/SuperMFE/xyt_1", "C:/Users/Alex/Downloads/SuperMFE/SuperMFE/grid_x.txt", "C:/Users/Alex/Downloads/SuperMFE/SuperMFE/grid_y.txt", "C:/Users/Alex/Downloads/SuperMFE/SuperMFE/grid_t.txt", condition1);

    system("pause");
    return 0;
}