Trud

#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <semaphore.h>
#include <sys/time.h>
#include <cmath>
#include <iostream>

#define DX 0.5
#define DY 0.5

#define DXR 0.25
#define DYR 0.25

#define BORDER_TEMPERATURE 200

#define HEIGHT 8
#define WIDTH 12
#define XPL 3
#define YPL 5

#define PSX ((WIDTH - XPL) / (2 * DX))
#define PEX ((WIDTH + XPL) / (2 * DX))
#define PSY ((HEIGHT - YPL) / (2 * DY))
#define PEY ((HEIGHT + YPL) / (2 * DY))
#define PSXR ((WIDTH - XPL) / (2 * DXR))
#define PEXR ((WIDTH + XPL) / (2 * DXR))
#define PSYR ((HEIGHT - YPL) / (2 * DYR))
#define PEYR ((HEIGHT + YPL) / (2 * DYR))

#define CALC_TIME 200
#define delta 4.0

using namespace std;

double fill_after(double ***T, int w, int h, int x, int y, double dx, double dy, int psx, int pex, int psy, int pey, int time, double DT)
{
    // для ГУ (левая правая стенка)
    if (x == 0 || x == w - 1)
    {
        return BORDER_TEMPERATURE;
    }
    // для верхней и нижней стенки ГУ 2-го рода
    else if (y == h - 1)
    {
        return T[time + 1][x][y - 1];
    }
    else if (y == 0)
    {
        return T[time][x][y + 1];
    }
    // если попадаем в отверстие
    else if (x < pex && x > psx && y < pey && y > psy)
    {
        return 0;
    }
    // левая стенка отверстия
    else if (x == psx && y < pey && y > psy)
    {
        return T[time][x - 1][y] / (sqrt(dx * dx + dy * dy) + 1);
    }
    // правая стенка отверстия
    else if (x == pex && y < pey && y > psy)
    {
        return T[time][x + 1][y] / (sqrt(dx * dx + dy * dy) + 1);
    }
    // нижняя стенка отверстия
    else if (x < pex && x > psx && y == psy)
    {
        return T[time][x][y - 1] / (sqrt(dx * dx + dy * dy) + 1);
    }
    // верхняя сторона отверстия
    else if (x < pex && x > psx && y == pey)
    {
        return T[time][x][y + 1] / (sqrt(dx * dx + dy * dy) + 1);
    }
    // все остальные узлы по формуле трехмерной явной разностной схеме
    else
    {
        return DT * ((T[time][x + 1][y] - 2 * T[time][x][y] + T[time][x - 1][y]) / (dx * dx) + (T[time][x][y + 1] - 2 * T[time][x][y] + T[time][x][y - 1]) / (dy * dy)) + T[time][x][y];
    }
}

int main()
{
    // для шага 0.5
    int h = int(HEIGHT / DY) + 1;
    int w = int(WIDTH / DX) + 1;

    // для шага 0.25
    int hr = int(HEIGHT / DYR) + 1;
    int wr = int(WIDTH / DXR) + 1;

    double maxder = 0;
    double **MP;

    // матрица для работы с шагом 0.5
    double **P;

    P = new double *[w];

    // начальное заполнение матриц
    for (int x = 0; x < w; x++)
    {
        P[x] = new double[h];
        for (int y = 0; y < h; y++)
        {
            P[x][y] = 0;
            if (x == 0 || x == w - 1)
            {
                P[x][y] = BORDER_TEMPERATURE;
            }
        }
    }

    MP = new double *[w - 1];

    for (int x = 0; x < w - 1; x++)
    {
        MP[x] = new double[h - 1];
        for (int y = 0; y < h - 1; y++)
        {
            MP[x][y] = 0;
        }
    }

    for (int x = 1; x < w - 1; x++)
    {
        for (int y = 1; y < h - 1; y++)
        {
            if (y == h - 2)
                MP[x - 1][y] = ((P[x + 1][y] - 2 * P[x][y] + P[x - 1][y]) / (DX * DX) + (P[x][y + 1] - 2 * P[x][y] + P[x][y - 1]) / (DY * DY));
            else if (x == w - 2)
                MP[x][y - 1] = ((P[x + 1][y] - 2 * P[x][y] + P[x - 1][y]) / (DX * DX) + (P[x][y + 1] - 2 * P[x][y] + P[x][y - 1]) / (DY * DY));
            else
                MP[x - 1][y - 1] = ((P[x + 1][y] - 2 * P[x][y] + P[x - 1][y]) / (DX * DX) + (P[x][y + 1] - 2 * P[x][y] + P[x][y - 1]) / (DY * DY));
        }
    }

    for (int x = 0; x < h - 1; x++)
    {
        for (int y = 0; y < w - 1; y++)
        {
            if (MP[x][y] > maxder)
                maxder = MP[x][y];
        }
    }



    double DT = (double)(delta / maxder);
    int MAX_TIME = (int)(CALC_TIME / DT);
    double ***T = new double **[MAX_TIME];
    double ***TR = new double **[MAX_TIME];

    // создание и инициализация трехмерной матрицы
    for (int i = 0; i < MAX_TIME; i++)
    {
        T[i] = new double *[w];
        for (int j = 0; j < w; j++)
        {
            T[i][j] = new double[h];
        }
    }

    for (int x = 0; x < w; x++)
    {
        for (int y = 0; y < h; y++)
        {
            T[0][x][y] = 0;
            if (x == 0 || x == w - 1)
            {
                T[0][x][y] = BORDER_TEMPERATURE;
            }
        }
    }

    for (int i = 0; i < MAX_TIME; i++)
    {
        TR[i] = new double *[wr];
        for (int j = 0; j < wr; j++)
        {
            TR[i][j] = new double[hr];
        }
    }

    for (int x = 0; x < wr; x++)
    {
        for (int y = 0; y < hr; y++)
        {
            TR[0][x][y] = 0;
            if (x == 0 || x == wr - 1)
            {
                TR[0][x][y] = BORDER_TEMPERATURE;
            }
        }
    }

    int timeCount = MAX_TIME - 1;

    for (int time = 0; time < timeCount; time++)
    {
        for (int x = 0; x < w; x++)
        {
            for (int y = 0; y < h; y++)
            {
                T[time + 1][x][y] = fill_after(T, w, h, x, y, DX, DY, PSX, PEX, PSY, PEY, time, DT);
            }
        }
    }

    cout << "h" << endl;

    FILE *out;
    out = fopen("out1.txt", "w");

    if (out == NULL)
    {
        perror("Ошибка открытия файла 1\n");
        exit(EXIT_FAILURE);
    }

    for (int y = h - 1; y >= 0; y--)
    {
        for (int x = 0; x < w; x++)
        {
            printf("%-6.2lf|", T[MAX_TIME - 1][x][y]);
            fprintf(out, "%d %d %lf\n", x, y, T[MAX_TIME - 1][x][y]);
        }
        printf("\n");
        fprintf(out, "\n");
    }

    for (int time = 0; time < timeCount; time++)
    {
        for (int x = 0; x < wr; x++)
        {
            for (int y = 0; y < hr; y++)
            {
                TR[time + 1][x][y] = fill_after(TR, wr, hr, x, y, DXR, DYR, PSXR, PEXR, PSYR, PEYR, time, DT);
            }
        }
    }

    printf("\n");

    FILE *out2;
    out2 = fopen("out2.txt", "w");

    if (out2 == NULL)
    {
        perror("Ошибка открытия файла 2\n");
        exit(EXIT_FAILURE);
    }

    for (int y = hr - 1; y >= 0; y--)
    {
        for (int x = 0; x < wr; x++)
        {
            printf("%-6.2lf|", TR[MAX_TIME - 1][x][y]);
            fprintf(out2, "%d %d %lf\n", x, y, TR[MAX_TIME - 1][x][y]);
        }
        printf("\n");
        fprintf(out2, "\n");
    }

    for (int x = 0; x < w; x++)
    {
        for (int y = 0; y < h; y++)
        {
            if (x == 0 || x == w - 1)
                P[x][y] = BORDER_TEMPERATURE;
            else if (y == h - 1)
                P[x][y] = P[x][y - 1];
            else if (y == 1)
                P[x][y] = P[x][y - 1];
            else
                P[x][y] = ((4 * T[MAX_TIME - 1][x][y]) - (TR[MAX_TIME - 1][2 * x][2 * y])) / 3;
        }
    }

    printf("\n");
    FILE *outr;
    outr = fopen("outr.txt", "w");

    if (outr == NULL)
    {
        perror("Ошибка открытия файла 3\n");
        exit(EXIT_FAILURE);
    }

    for (int y = h - 1; y >= 0; y--)
    {
        for (int x = 0; x < w; x++)
        {
            printf("%-6.2lf|", P[x][y]);
            fprintf(outr, "%d %d %lf\n", x, y, P[x][y]);
        }
        printf("\n");
        fprintf(outr, "\n");
    }
    return 0;
}