API tools faq
paste
Gistrec

МКЭ Курсовая

Gistrec
Dec 3rd, 2019
327
0
Never
1
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 31.15 KB | None | 0 0
raw download clone embed print report
  1. #include <iostream>
  2. #include <fstream>
  3. #include <functional>
  4. #include <vector>
  5. #include <string>
  6. #include <tuple>
  7.  
  8. using std::cout;
  9. using std::ifstream;
  10. using std::ofstream;
  11. using std::endl;
  12. using std::function;
  13. using std::vector;
  14. using std::string;
  15. using std::ostream;
  16. using std::istream;
  17. using std::swap;
  18. using std::tuple;
  19. using std::make_tuple;
  20. using std::tie;
  21.  
  22. #define DIM_LOCAL 4
  23.  
  24. //--------------------------------------------------------------------------------
  25. // Структура для хранение тестировочных данных
  26. // U        - истиннове решение (x, y, t)
  27. // f        - функция правой части (x, y, t)
  28. // sigma    - сигма из уравнения
  29. // hee      - хи из уравнения
  30. // lambda   - лямбда из уравнения
  31. //--------------------------------------------------------------------------------
  32. struct TestData {
  33.     function<double(double, double, double)> u;
  34.     function<double(double, double, double)> f;
  35.     double sigma;
  36.     double hee;
  37.     double lambda;
  38. };
  39.  
  40. //--------------------------------------------------------------------------------
  41. // Структура для хранения простейших первых граничных условий,
  42. // прямоугольная область, на пересечениях границ функции должны быть равны
  43. // leftXCondition   - левая граница области(x = x0)
  44. // rightXCondition  - правая граница области(x = x_max)
  45. // leftYCondition   - нижняя граница области(y = y0)
  46. // rightYCondition  - верхняя граница области(y = y_max)
  47. //--------------------------------------------------------------------------------
  48. struct Temp1stBoundCondition {
  49.     function<double(double, double, double)> leftXCondition;
  50.     function<double(double, double, double)> rightXCondition;
  51.     function<double(double, double, double)> leftYCondition;
  52.     function<double(double, double, double)> rightYCondition;
  53. };
  54.  
  55. //--------------------------------------------------------------------------------
  56. // Функция загрузки и конструирования сетки
  57. // Формат данных в потоке: начальная точка, конечная точка,
  58. //  количество точек с начальной и конечной, коэффициент растяжения/сжатия
  59. // in       - входной поток
  60. // grid     - выход, координаты узлов сетки
  61. //--------------------------------------------------------------------------------
  62. void uploadGrid(istream & in, vector<double> & grid) {
  63.     double start, end, n, k;
  64.     double h;
  65.     in >> start >> end >> n >> k;
  66.     n = n - 1;
  67.     if (!grid.empty()) grid.clear();
  68.  
  69.     grid.push_back(start);
  70.     if (k != 1.) {
  71.         h = (end - start) * (1 - k) / (1 - pow(k, n));
  72.         for (int i = 1; i < n; i++) {
  73.             grid.push_back(grid.back() + h);
  74.             h *= k;
  75.         }
  76.     }
  77.     else {
  78.         h = (end - start) / (n);
  79.         for (int i = 1; i < n; i++) {
  80.             grid.push_back(grid.back() + h);
  81.         }
  82.     }
  83.     grid.push_back(end);
  84. }
  85.  
  86. //--------------------------------------------------------------------------------
  87. // Постройка сеток для двумерной нестационарной задачи
  88. // xGrid        - выход, сетка по иксу
  89. // yGrid        - выход, сетка по игреку
  90. // tGrid        - выход, сетка по времени
  91. // xGridName    - имя файла с параметрами сетки по иксу
  92. // yGridName    - имя файла с параметрами сетки по игреку
  93. // tGridName    - имя файла с параметрами сетки по времени
  94. //--------------------------------------------------------------------------------
  95. void constructGrids(vector<double> & xGrid, vector<double> & yGrid, vector<double> & tGrid,
  96.                     const string & xGridName = "grid_x.txt", const string & yGridName = "grid_y.txt", const string & tGridName = "grid_t.txt") {
  97.     ifstream in(xGridName);
  98.     uploadGrid(in, xGrid);
  99.     in.close();
  100.     in.open(yGridName);
  101.     uploadGrid(in, yGrid);
  102.     in.close();
  103.     in.open(tGridName);
  104.     uploadGrid(in, tGrid);
  105.     in.close();
  106. }
  107.  
  108. //--------------------------------------------------------------------------------
  109. // Ядро локальной матрицы Жёсткости(с множителем l/6*hy/hx)
  110. //--------------------------------------------------------------------------------
  111. double kernGLocalYX[DIM_LOCAL][DIM_LOCAL] = {   { 2, -2, 1, -1},
  112.                                                 { -2, 2, -1, 1},
  113.                                                 { 1, -1, 2, -2},
  114.                                                 { -1, 1, -2, 2} };
  115.  
  116. //--------------------------------------------------------------------------------
  117. // Ядро локальной матрицы Жёсткости(с множителем l/6*hx/hy)
  118. //--------------------------------------------------------------------------------
  119. double kernGLocalXY[DIM_LOCAL][DIM_LOCAL] = {   { 2, 1, -2, -1 },
  120.                                                 { 1, 2, -1, -2 },
  121.                                                 { -2, -1, 2, 1 },
  122.                                                 { -1, -2, 1, 2 } };
  123.  
  124. //--------------------------------------------------------------------------------
  125. // Ядро локальной матрицы Масс(Ядро локальной матрицы C)
  126. //--------------------------------------------------------------------------------
  127. double kernMLocal[DIM_LOCAL][DIM_LOCAL] = { { 4, 2, 2, 1 },
  128.                                             { 2, 4, 1, 2 },
  129.                                             { 2, 1, 4, 2 },
  130.                                             { 1, 2, 2, 4 } };
  131.  
  132. //--------------------------------------------------------------------------------
  133. // Получение компоненты локального вектора b
  134. // index    - элемент вектора b(счёт с нуля)
  135. // f        - функция правой части (x, y, t)
  136. // x1       - начальная координата по иксу конечного элемента
  137. // x2       - конечная координата по иксу конечного элемента
  138. // y1       - начальная координата по игреку конечного элемента
  139. // y2       - конечная координата по игреку конечного элемента
  140. // t        - значение времени текущего временного слоя
  141. //--------------------------------------------------------------------------------
  142. double getBLocal(const size_t index, const function<double(double, double, double)> & f,
  143.                 const double x1, const double x2, const double y1, const double y2, const double t) {
  144.     double result = 0;
  145.  
  146.     result += kernMLocal[index][0] * f(x1, y1, t);
  147.     result += kernMLocal[index][1] * f(x2, y1, t);
  148.     result += kernMLocal[index][2] * f(x1, y2, t);
  149.     result += kernMLocal[index][3] * f(x2, y2, t);
  150.  
  151.     return result;// *(x2 - x1) * (y2 - y1) / 36 * (x2 - x1) * (y2 - y1) / 36;
  152. }
  153.  
  154. //--------------------------------------------------------------------------------
  155. // Вставляет локальную матрицу в глобальную
  156. // local    - локальная матрица
  157. // global   - глобальная матрица
  158. // i1       - первый глобальный индекс вставки
  159. // i2       - второй глобальный индекс вставки
  160. // i3       - третий глобальный индекс вставки
  161. // i4       - четвёртый глобальный индекс вставки
  162. //--------------------------------------------------------------------------------
  163. void insertLocalToGlobal(const double (*local)[DIM_LOCAL], double ** global, const double mult,
  164.     const size_t i1, const size_t i2, const size_t i3, const size_t i4) {
  165.  
  166.     vector<size_t> indexes = { i1, i2, i3, i4 };
  167.     for (int i = 0; i < DIM_LOCAL; i++) {
  168.         for (int j = 0; j < DIM_LOCAL; j++) {
  169.             global[indexes[i]][indexes[j]] += mult * local[i][j];
  170.         }
  171.     }
  172. }
  173.  
  174. //--------------------------------------------------------------------------------
  175. // Вставляет локальную матрицу в глобальную
  176. // local    - локальная матрица
  177. // global   - глобальная матрица
  178. // i1       - первый глобальный индекс вставки
  179. // i2       - второй глобальный индекс вставки
  180. // i3       - третий глобальный индекс вставки
  181. // i4       - четвёртый глобальный индекс вставки
  182. //--------------------------------------------------------------------------------
  183. void insertLocalToGlobal(const double * local, double * global, const double mult,
  184.     const size_t i1, const size_t i2, const size_t i3, const size_t i4) {
  185.  
  186.     global[i1] += mult * local[0];
  187.     global[i2] += mult * local[1];
  188.     global[i3] += mult * local[2];
  189.     global[i4] += mult * local[3];
  190. }
  191.  
  192. //--------------------------------------------------------------------------------
  193. // Вспомогательная функция для вычисление глобального индекса первого узла по
  194. // номеру конечного элемента
  195. // numberFE - номер конечного элемента(счёт с 0)
  196. // sizeX    - количество узлов по иксу
  197. //--------------------------------------------------------------------------------
  198. inline size_t globalStartIndex(const size_t numberFE, const size_t sizeX) {
  199.     return numberFE + static_cast<size_t>(numberFE / (sizeX - 1));
  200. }
  201.  
  202. //--------------------------------------------------------------------------------
  203. // Вспомогательная функция для вычисление индекса икса по глобальному индексу узла
  204. // numberFE - номер конечного элемента(счёт с 0)
  205. // sizeX    - количество узлов по иксу
  206. //--------------------------------------------------------------------------------
  207. inline size_t indexXbyGlobalIndex(const size_t numberFE, const size_t sizeX) {
  208.     return globalStartIndex(numberFE,sizeX) % sizeX;
  209. }
  210.  
  211. //--------------------------------------------------------------------------------
  212. // Вспомогательная функция для вычисление индекса игрека по глобальному индексу
  213. // узла
  214. // numberFE - номер конечного элемента(счёт с 0)
  215. // sizeX    - количество узлов по иксу
  216. //--------------------------------------------------------------------------------
  217. inline size_t indexYbyGlobalIndex(const size_t numberFE, const size_t sizeX) {
  218.     return static_cast<size_t>(globalStartIndex(numberFE, sizeX) / sizeX);
  219. }
  220.  
  221. //--------------------------------------------------------------------------------
  222. // Сборка глобальной матрицы Масс(гамма)
  223. // globalM  - глобальная матрица Масс
  224. // xGrid    - сетка по иксу
  225. // yGrid    - сетка по игреку
  226. // gamma    - гамма(множитель)
  227. //--------------------------------------------------------------------------------
  228. void constructGlobalM(double ** globalM, const vector<double> & xGrid,
  229.     const vector<double> & yGrid, const double gamma) {
  230.  
  231.     size_t countFE = (xGrid.size() - 1) * (yGrid.size() - 1);
  232.     double dx, dy;
  233.     double xStart, xEnd;
  234.     double yStart, yEnd;
  235.     double globalStart;
  236.  
  237.     for (int elemNum = 0; elemNum < countFE; elemNum++) {
  238.         globalStart = globalStartIndex(elemNum, xGrid.size());
  239.         xStart = xGrid[indexXbyGlobalIndex(elemNum, xGrid.size())];
  240.         xEnd = xGrid[indexXbyGlobalIndex(elemNum, xGrid.size()) + 1];
  241.         yStart = yGrid[indexYbyGlobalIndex(elemNum, xGrid.size())];
  242.         yEnd = yGrid[indexYbyGlobalIndex(elemNum, xGrid.size()) + 1];
  243.  
  244.         dx = xEnd - xStart;
  245.         dy = yEnd - yStart;
  246.  
  247.         insertLocalToGlobal(kernMLocal, globalM, gamma * dx * dy / 36,
  248.             globalStart, globalStart + 1, globalStart + xGrid.size(), globalStart + xGrid.size() + 1);
  249.     }
  250. }
  251.  
  252. //--------------------------------------------------------------------------------
  253. // Сборка глобальный матрицы Жёсткости(лямбда)
  254. // globalM  - глобальная матрица Масс
  255. // xGrid    - сетка по иксу
  256. // yGrid    - сетка по игреку
  257. // lambda   - лямбда(множитель)
  258. //--------------------------------------------------------------------------------
  259. void constructGlobalG(double ** globalG, const vector<double> & xGrid,
  260.     const vector<double> & yGrid, const double lambda) {
  261.  
  262.     size_t countFE = (xGrid.size() - 1) * (yGrid.size() - 1);
  263.     double dx, dy;
  264.     double xStart, xEnd;
  265.     double yStart, yEnd;
  266.     double globalStart;
  267.  
  268.     for (int elemNum = 0; elemNum < countFE; elemNum++) {
  269.         globalStart = globalStartIndex(elemNum, xGrid.size());
  270.         xStart = xGrid[indexXbyGlobalIndex(elemNum, xGrid.size())];
  271.         xEnd = xGrid[indexXbyGlobalIndex(elemNum, xGrid.size()) + 1];
  272.         yStart = yGrid[indexYbyGlobalIndex(elemNum, xGrid.size())];
  273.         yEnd = yGrid[indexYbyGlobalIndex(elemNum, xGrid.size()) + 1];
  274.  
  275.         dx = xEnd - xStart;
  276.         dy = yEnd - yStart;
  277.    
  278.         insertLocalToGlobal(kernGLocalYX, globalG, lambda * dy / 6. / dx,
  279.             globalStart, globalStart + 1, globalStart + xGrid.size(), globalStart + xGrid.size() + 1);
  280.         insertLocalToGlobal(kernGLocalXY, globalG, lambda * dx / 6. / dy,
  281.             globalStart, globalStart + 1, globalStart + xGrid.size(), globalStart + xGrid.size() + 1);
  282.     }
  283. }
  284.  
  285. //--------------------------------------------------------------------------------
  286. // Сборка глобального вектора правой части
  287. // globalB  - глобальный вектор B
  288. // xGrid    - сетка по иксу
  289. // yGrid    - сетка по игреку
  290. // f        - функция правой части
  291. // t        - текущий временной слой
  292. //--------------------------------------------------------------------------------
  293. void constructGlobalB(double * globalB, const vector<double> & xGrid,
  294.     const vector<double> & yGrid, const function<double(double,double,double)> f,
  295.     const double t) {
  296.  
  297.     size_t countFE = (xGrid.size() - 1) * (yGrid.size() - 1);
  298.     double dx, dy;
  299.     double xStart, xEnd;
  300.     double yStart, yEnd;
  301.     double globalStart;
  302.  
  303.     double localRightPart[DIM_LOCAL];
  304.  
  305.     for (int elemNum = 0; elemNum < countFE; elemNum++) {
  306.         globalStart = globalStartIndex(elemNum, xGrid.size());
  307.         xStart = xGrid[indexXbyGlobalIndex(elemNum, xGrid.size())];
  308.         xEnd = xGrid[indexXbyGlobalIndex(elemNum, xGrid.size()) + 1];
  309.         yStart = yGrid[indexYbyGlobalIndex(elemNum, xGrid.size())];
  310.         yEnd = yGrid[indexYbyGlobalIndex(elemNum, xGrid.size()) + 1];
  311.  
  312.         dx = xEnd - xStart;
  313.         dy = yEnd - yStart;
  314.  
  315.         for (int i = 0; i < DIM_LOCAL; i++)
  316.             localRightPart[i] = getBLocal(i, f, xStart, xEnd, yStart, yEnd, t);
  317.  
  318.         insertLocalToGlobal(localRightPart, globalB, dx * dy / 36.,
  319.             globalStart, globalStart + 1, globalStart + xGrid.size(), globalStart + xGrid.size() + 1);
  320.     }
  321. }
  322.  
  323. //--------------------------------------------------------------------------------
  324. // Решатель СЛАУ методом Гаусс.
  325. // globalA  - глобальная матрица A
  326. // globalB  - вектор правой части(вход), результат решения(выход)
  327. // n        - размернойсть квадратной матрицы и вектора правой части
  328. //--------------------------------------------------------------------------------
  329. void solveSLAE(const double * const * globalA, double *globalB, size_t n) {
  330.     double tmp;
  331.     double **A = new double *[n];
  332.     for (int i = 0; i < n; i++) {
  333.         A[i] = new double[n];
  334.         for (int j = 0; j < n; j++)
  335.             A[i][j] = globalA[i][j];
  336.     }
  337.  
  338.     for (int i = 0; i < n; i++) {
  339.         tmp = A[i][i];
  340.         for (int j = i; j < n; j++)
  341.             A[i][j] /= tmp;
  342.         globalB[i] /= tmp;
  343.         for (int j = i + 1; j < n; j++) {
  344.             tmp = A[j][i];
  345.             for (int k = i; k < n; k++)
  346.                 A[j][k] -= tmp * A[i][k];
  347.             globalB[j] -= tmp * globalB[i];
  348.         }
  349.         for (int j = i - 1; j >= 0; j--) {
  350.             tmp = A[j][i];
  351.             for (int k = i; k < n; k++)
  352.                 A[j][k] -= tmp * A[i][k];
  353.             globalB[j] -= tmp * globalB[i];
  354.         }
  355.     }
  356.  
  357.     for (int i = 0; i < n; i++)
  358.         delete[] A[i];
  359.     delete[] A;
  360. }
  361.  
  362. //--------------------------------------------------------------------------------
  363. // Применяем первый краевые условия
  364. // globalA         - глобальная матрица для СЛАУ
  365. // globalRightPart - глобальный вектор правой части для СЛАУ
  366. // xGrid           - сетка по иксу
  367. // yGrid           - сетка по игреку
  368. // condition       - структура с краевыми условиями первого рода
  369. // t               - текущий временной слой
  370. //--------------------------------------------------------------------------------
  371. void apply1stCondition(double **globalA, double *globalRightPart,
  372.     const vector<double> & xGrid, const vector<double> & yGrid,
  373.     const Temp1stBoundCondition & condition, const double t) {
  374.  
  375.     size_t countNode = (xGrid.size()) * (yGrid.size());
  376.  
  377.     // Нижняя граница. y = y_min, x = [x_min, x_max]
  378.     for (int i = 0; i < xGrid.size(); i++) {
  379.         for (int j = 0; j < countNode; j++) {
  380.             globalA[i][j] = 0;
  381.         }
  382.         globalA[i][i] = 1.;
  383.  
  384.         globalRightPart[i] = condition.leftYCondition(xGrid[i], yGrid[0], t);
  385.     }
  386.     // Верхняя граница. y = y_max, x = [x_min, x_max]
  387.     for (int xIndex = 0, curElem = countNode - xGrid.size(); curElem < countNode; xIndex++, curElem++) {
  388.         for (int j = 0; j < countNode; j++) {
  389.             globalA[curElem][j] = 0;
  390.         }
  391.         globalA[curElem][curElem] = 1.;
  392.  
  393.         globalRightPart[curElem] = condition.rightYCondition(xGrid[xIndex], yGrid.back(), t);
  394.     }
  395.     // Левая граница. x = x_min, y = [y_min, y_max]
  396.     for (int yIndex = 0, curElem = 0; curElem < countNode; yIndex++, curElem += xGrid.size()) {
  397.         for (int j = 0; j < countNode; j++) {
  398.             globalA[curElem][j] = 0;
  399.         }
  400.         globalA[curElem][curElem] = 1.;
  401.  
  402.         globalRightPart[curElem] = condition.leftXCondition(xGrid[0], yGrid[yIndex], t);
  403.     }
  404.     // Правая граница. x = x_max, y = [y_min, y_max]
  405.     for (int yIndex = 0, curElem = xGrid.size() - 1; curElem < countNode; yIndex++, curElem += xGrid.size()) {
  406.         for (int j = 0; j < countNode; j++) {
  407.             globalA[curElem][j] = 0;
  408.         }
  409.         globalA[curElem][curElem] = 1.;
  410.  
  411.         globalRightPart[curElem] = condition.rightXCondition(xGrid.back(), yGrid[yIndex], t);
  412.     }
  413. }
  414.  
  415. //--------------------------------------------------------------------------------
  416. // Постройка весового вектора по истинному решению
  417. // qi       - глобальный вектор весов
  418. // ti       - значение времени на временном слое
  419. // xGrid    - сетка по иксу
  420. // yGrid    - сетка по игреку
  421. // u        - истинное решение
  422. //--------------------------------------------------------------------------------
  423. void constructQiSyntetic(double * qi, const double ti, const vector<double> & xGrid,
  424.     const vector<double> & yGrid, const function<double(double, double, double)> & u) {
  425.  
  426.     for (int curY = 0, i = 0; curY < yGrid.size(); curY++) {
  427.         for (int curX = 0; curX < xGrid.size(); curX++, i++) {
  428.             qi[i] = u(xGrid[curX], yGrid[curY], ti);
  429.         }
  430.     }
  431. }
  432.  
  433. //--------------------------------------------------------------------------------
  434. // Сборка начальных слоёв для трёхточечной схемы
  435. // q_2      - выход, вектор весов соответствующий t = t_0
  436. // q_1      - выход, вектор весов соответствующий t = t_1
  437. // t_2      - значение времени t_0
  438. // t_1      - значение времени t_1
  439. // xGrid    - сетка по иксу
  440. // yGrid    - сетка по игреку
  441. // u        - истинное решение
  442. //--------------------------------------------------------------------------------
  443. void constructStartCondition(double * q_2, double * q_1, const double t_2, const double t_1,
  444.     const vector<double> & xGrid, const vector<double> & yGrid,
  445.     const function<double(double, double, double)> & u) {
  446.  
  447.     constructQiSyntetic(q_2, t_2, xGrid, yGrid, u);
  448.     constructQiSyntetic(q_1, t_1, xGrid, yGrid, u);
  449.  
  450. }
  451.  
  452. //--------------------------------------------------------------------------------
  453. // Правило сборки глобальное матрицы СЛАУ для гиперболического уравнения со схемой
  454. // Кранка-Николсона(трёхточечная неявная по времени)
  455. // globalG      - глобальная марица жёсткости уже умноженная на лямбду
  456. // globalMhi    - глобальная матрица масс умноженная на хи
  457. // globalMsig   - глобальная матрица масс умноженная на сигму
  458. // globalA      - глобальная матрица СЛАУ
  459. // dt           - шаг по времени
  460. // dim          - количество узлов в глобальной сетки(размерность весового вектора)
  461. //--------------------------------------------------------------------------------
  462. void constructGlobalA(const double * const * globalG, const double * const * globalMhi,
  463.     const double * const * globalMsig, double * const * globalA, const double dt, const size_t dim) {
  464.  
  465.     for (int row = 0; row < dim; row++) {
  466.         for (int col = 0; col < dim; col++) {
  467.             globalA[row][col] = 1. / dt / dt * globalMhi[row][col]
  468.                 + 1. / 2. / dt * globalMsig[row][col] + 1. / 2. * globalG[row][col];
  469.         }
  470.     }
  471. }
  472.  
  473. //--------------------------------------------------------------------------------
  474. // TODO: и это дозаполнять
  475. //--------------------------------------------------------------------------------
  476. void constructGlobalRightPart(const double * globalB, const double * globalB_2,
  477.     const double * const * globalMhi, const double * const * globalMsig, const double *const * globalG,
  478.     const double * qj_1, const double * qj_2, const double dt, double * globalRight,
  479.     const size_t dim) {
  480.    
  481.     double temp;
  482.  
  483.     for (int curLine = 0; curLine < dim; curLine++) {
  484.         temp = 0;
  485.         // + 1/2 bj
  486.         globalRight[curLine] = 1. / 2. * globalB[curLine];
  487.         // + 1/2 bj_2
  488.         globalRight[curLine] += 1. / 2. * globalB_2[curLine];
  489.         // + 2. / dt^2 * Mhi * qj-1
  490.         for (int sumIndex = 0; sumIndex < dim; sumIndex++) {
  491.             temp += globalMhi[curLine][sumIndex] * qj_1[sumIndex];
  492.         }
  493.         globalRight[curLine] += 2. / dt / dt * temp;
  494.         temp = 0;
  495.         // - 1./dt^2 * Mhi * qj-2
  496.         for (int sumIndex = 0; sumIndex < dim; sumIndex++) {
  497.             temp += globalMhi[curLine][sumIndex] * qj_2[sumIndex];
  498.         }
  499.         globalRight[curLine] -= 1. / dt / dt * temp;
  500.         temp = 0;
  501.         // + 1/2dt * Msig * qj-2
  502.         for (int sumIndex = 0; sumIndex < dim; sumIndex++) {
  503.             temp += globalMsig[curLine][sumIndex] * qj_2[sumIndex];
  504.         }
  505.         globalRight[curLine] += 1. / 2. / dt * temp;
  506.         temp = 0;
  507.         // - 1/2 * G qj-2
  508.         for (int sumIndex = 0; sumIndex < dim; sumIndex++) {
  509.             temp += globalG[curLine][sumIndex] * qj_2[sumIndex];
  510.         }
  511.         globalRight[curLine] -= 1. / 2. * temp;
  512.     }
  513.  
  514. }
  515.  
  516.  
  517. //--------------------------------------------------------------------------------
  518. // Вывод в файл правой части
  519. // filename     - имя файла
  520. // rightPart    - вектор правой части
  521. // dim          - размерность вектора
  522. //--------------------------------------------------------------------------------
  523. void outputRightPart(const string & fileName, const double * rightPart, const size_t dim) {
  524.     ofstream out(fileName);
  525.     for (int i = 0; i < dim; i++)
  526.         out << rightPart[i] << endl;
  527.     out.close();
  528. }
  529.  
  530. //--------------------------------------------------------------------------------
  531. // TODO: Заполнить
  532. //--------------------------------------------------------------------------------
  533. void calcErr(const string & fileTruth, const string & curFile, const string & outFile) {
  534.     ifstream inTruth(fileTruth);
  535.     ifstream inCur(curFile);
  536.     ofstream out(outFile);
  537.  
  538.     while (!inTruth.eof()) {
  539.         out << inTruth.get() - inCur.get() << endl;
  540.     }
  541.  
  542.     out.close();
  543.     inCur.close();
  544.     inTruth.close();
  545. }
  546.  
  547. //--------------------------------------------------------------------------------
  548. // TODO: Заполнить
  549. //--------------------------------------------------------------------------------
  550. void swapVectors(double * v1, double * v2, const size_t dim) {
  551.     for (int i = 0; i < dim; i++)
  552.         swap(v1[i], v2[i]);
  553. }
  554.  
  555. //--------------------------------------------------------------------------------
  556. // TODO: Заполнить и это не забыть
  557. //--------------------------------------------------------------------------------
  558. void clearVector(double * v, const size_t dim) {
  559.     for (int i = 0; i < dim; i++)
  560.         v[i] = 0;
  561. }
  562.  
  563. //--------------------------------------------------------------------------------
  564. // TODO: Заполнить и это не забыть
  565. //--------------------------------------------------------------------------------
  566. void testing(const TestData & test, const string & testName, const string & xFileName,
  567.              const string & yFileName, const string & tFileName,
  568.              const Temp1stBoundCondition & condition) {
  569.  
  570.     vector<double> xGrid, yGrid, tGrid;
  571.     constructGrids(xGrid, yGrid, tGrid, xFileName, yFileName, tFileName);
  572.  
  573.     size_t countFE = (xGrid.size() - 1) * (yGrid.size() - 1);
  574.     size_t countNode = xGrid.size() * yGrid.size();
  575.  
  576.     double ** globalMhi = new double*[countNode];
  577.     for (int i = 0; i < countNode; i++)
  578.         globalMhi[i] = new double[countNode]();
  579.  
  580.     double ** globalMsig = new double*[countNode];
  581.     for (int i = 0; i < countNode; i++)
  582.         globalMsig[i] = new double[countNode]();
  583.  
  584.     double ** globalG = new double*[countNode];
  585.     for (int i = 0; i < countNode; i++)
  586.         globalG[i] = new double[countNode]();
  587.  
  588.     double ** globalA = new double*[countNode];
  589.     for (int i = 0; i < countNode; i++)
  590.         globalA[i] = new double[countNode]();
  591.  
  592.     double * curGlobalB = new double[countNode];
  593.     double * curGlobalB_2 = new double[countNode];
  594.     double * curQ = new double[countNode];
  595.     double * q_1 = new double[countNode];
  596.     double * q_2 = new double[countNode];
  597.     double * globalRight = new double[countNode];
  598.  
  599.  
  600.     constructStartCondition(q_2, q_1, tGrid[0], tGrid[1], xGrid, yGrid, test.u);
  601.     constructGlobalM(globalMhi, xGrid, yGrid, test.hee);
  602.     constructGlobalM(globalMsig, xGrid, yGrid, test.sigma);
  603.     constructGlobalG(globalG, xGrid, yGrid, test.lambda);
  604.  
  605.     outputRightPart(testName + "\\t_" + std::to_string(tGrid[0]) + ".txt", q_2, countNode);
  606.     outputRightPart(testName + "\\t_" + std::to_string(tGrid[1]) + ".txt", q_1, countNode);
  607.  
  608.     double dt;
  609.     for (int curT = 2; curT < tGrid.size(); curT++) {
  610.         dt = tGrid[curT] - tGrid[curT - 1];
  611.         constructGlobalA(globalG, globalMhi, globalMsig, globalA, dt, countNode);
  612.         clearVector(curGlobalB, countNode);
  613.         clearVector(curGlobalB_2, countNode);
  614.  
  615.         constructGlobalB(curGlobalB, xGrid, yGrid, test.f, tGrid[curT]);
  616.         constructGlobalB(curGlobalB_2, xGrid, yGrid, test.f, tGrid[curT - 2]);
  617.         constructGlobalRightPart(curGlobalB, curGlobalB_2, globalMhi, globalMsig, globalG, q_1, q_2, dt, curQ, countNode);
  618.         apply1stCondition(globalA, curQ, xGrid, yGrid, condition, tGrid[curT]);
  619.         solveSLAE(globalA, curQ, countNode);
  620.  
  621.         swapVectors(q_2, q_1, countNode);
  622.         swapVectors(q_1, curQ, countNode);
  623.         constructQiSyntetic(curQ, tGrid[curT], xGrid, yGrid, test.u);
  624.  
  625.         outputRightPart(testName + "\\t_" + std::to_string(tGrid[curT]) + ".txt", q_1, countNode);
  626.         outputRightPart(testName + "\\truth_t_" + std::to_string(tGrid[curT]) + ".txt", curQ, countNode);
  627.         for (int i = 0; i < countNode; i++) {
  628.             curQ[i] -= q_1[i];
  629.         }
  630.         outputRightPart(testName + "\\err_t_" + std::to_string(tGrid[curT]) + ".txt", curQ, countNode);
  631.         outputRightPart(testName + "\\last_error.txt", curQ, countNode);
  632.     }
  633.  
  634.     double sum = 0.0;
  635.     for (int i = 0; i < countNode; i++) {
  636.         sum += curQ[i];
  637.     }
  638.     sum = sum / countNode;
  639.     std::cout << "Погрешность: " << sum << std::endl;
  640.  
  641.     delete[] globalRight;
  642.     delete[] q_2;
  643.     delete[] q_1;
  644.     delete[] curQ;
  645.     delete[] curGlobalB_2;
  646.     delete[] curGlobalB;
  647.  
  648.     for (int i = 0; i < countNode; i++)
  649.         delete[] globalA[i];
  650.     delete[] globalA;
  651.  
  652.     for (int i = 0; i < countNode; i++)
  653.         delete[] globalG[i];
  654.     delete[] globalG;
  655.  
  656.     for (int i = 0; i < countNode; i++)
  657.         delete[] globalMsig[i];
  658.     delete[] globalMsig;
  659.  
  660.     for (int i = 0; i < countNode; i++)
  661.         delete[] globalMhi[i];
  662.     delete[] globalMhi;
  663.  
  664. }
  665.  
  666.  
  667. //--------------------------------------------------------------------------------
  668. //
  669. //--------------------------------------------------------------------------------
  670.  
  671. int main() {
  672.    
  673.     // TODO: запустить тесты, отладить.
  674.     TestData test1;
  675.     test1.hee = 1;
  676.     test1.lambda = 6;//1;
  677.     test1.sigma = 1;
  678.     test1.u = [](double x, double y, double t) -> double {
  679.         return x * y*(t - 1);
  680.     };
  681.     test1.f = [](double x, double y, double t) -> double {
  682.         return x * y;
  683.     };
  684.     Temp1stBoundCondition condition1;
  685.     condition1.leftXCondition = [](double x, double y, double t) {
  686.         return  x * y * (t - 1);
  687.     };
  688.     condition1.leftYCondition = [](double x, double y, double t) {
  689.         return  x * y * (t - 1);
  690.     };
  691.     condition1.rightXCondition = [](double x, double y, double t) {
  692.         return x * y * (t - 1);
  693.     };
  694.     condition1.rightYCondition = [](double x, double y, double t) {
  695.         return  x * y * (t - 1);
  696.     };
  697.     //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);
  698.  
  699.  
  700.  
  701.     test1.hee = 1;
  702.     test1.lambda = 3;
  703.     test1.sigma = 2;
  704.     test1.u = [](double x, double y, double t) -> double {
  705.         return x * y * t;
  706.     };
  707.     test1.f = [](double x, double y, double t) -> double {
  708.         return 2 * x * y;
  709.     };
  710.     condition1.leftXCondition = test1.u;
  711.     condition1.leftYCondition = test1.u;
  712.     condition1.rightXCondition = test1.u;
  713.     condition1.rightYCondition = test1.u;
  714.     //testing(test1, "C:/Users/Alex/Downloads/SuperMFE/SuperMFE/xyt", "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);
  715.  
  716.      
  717.     test1.u = [](double x, double y, double t) -> double {
  718.         return x * x * x * y * y * t;
  719.     };
  720.     test1.f = [](double x, double y, double t) -> double {
  721.         return 2 * x * x * x * y * y - 18 * x * y * y * t - 6 * x * x * x * t;
  722.     };
  723.     condition1.leftXCondition = test1.u;
  724.     condition1.leftYCondition = test1.u;
  725.     condition1.rightXCondition = test1.u;
  726.     condition1.rightYCondition = test1.u;
  727.     testing(test1, "C:/Users/Alex/Downloads/SuperMFE/SuperMFE/x3y2t", "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);
  728.  
  729.  
  730.     test1.u = [](double x, double y, double t) -> double {
  731.         return x * x * x * x * x * x * y * t;
  732.     };
  733.     test1.f = [](double x, double y, double t) -> double {
  734.         return 2 * x * x * x * x * x * x * y - 90 * t * x * x * x * x * y;
  735.     };
  736.  
  737.     condition1.leftXCondition = test1.u;
  738.     condition1.leftYCondition = test1.u;
  739.     condition1.rightXCondition = test1.u;
  740.     condition1.rightYCondition = test1.u;
  741.     //testing(test1, "C:/Users/Alex/Downloads/SuperMFE/SuperMFE/x5yt", "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);
  742.  
  743.  
  744.     system("pause");
  745.     return 0;
  746. }
Comments
  • User was banned
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!