[cpp-qt] dd

Viewer

  1. #include <cmath>
  2. #include <fstream>
  3. #include <iostream>
  4. #include <vector>
  5. #include <set>
  6. #include <algorithm>
  7. #include <functional>
  8. #include "iomanip"
  9. using namespace std;
  10.  
  11. //узел 
  12. struct Node
  13. {
  14.     double r, z; // координаты узла
  15.     int globalNumber; // глобальный номер узла 
  16. };
  17.  
  18. // расчетная область
  19. struct Range
  20. {
  21. public:
  22.     double diffusionCoefficient; //коэф. диффизии
  23.     double gamma; // гамма коэф.
  24.     int rangeNumber;
  25.  
  26.     double GetValueF(Node node, double t)
  27.     {
  28.         switch (rangeNumber)
  29.         {
  30.         case 0:
  31.             return 20;
  32.         case 1:
  33.             return  24 * node.z + 16;
  34.         case 2:
  35.             return  24 * node.r + 16 - ((18 * t) / node.r);
  36.         case 3:
  37.             return (-36) * t + 24 * node.z * node.z + 16;
  38.         case 4:
  39.             return (-72) * t + 24 * node.r * node.r + 16;
  40.         case 5:
  41.             return 40 * t;
  42.         case 6:
  43.             return t * (48 * node.z + 32);
  44.         case 7:
  45.             return t * (48 * node.r + 32) - ((18 * t * t) / node.r);
  46.         case 8:
  47.             return t * (48 * node.z * node.z + 32) - 12 * t * t;
  48.         case 9:
  49.             return t * (48 * node.r * node.r + 32) - 24 * t * t;
  50.         default:
  51.             break;
  52.         }
  53.     }
  54. };
  55.  
  56. struct Triangle
  57. {
  58.     Node node[3]; // локальные узлы
  59.     Range range; // расченая области
  60.     vector<vector<double>> matrixMass; //матрица масс
  61.     vector<vector<double>> matrixStiffness; // матрица жесткости
  62.     vector<vector<double>> localMatrix; // локальная матрица 
  63.     vector<double> localVector; // локальный вектор 
  64.  
  65.     Triangle()
  66.     {
  67.         vector<vector<double>> matrixMass(3, vector<double>(3)); //матрица масс
  68.         vector<vector<double>> matrixStiffness(3, vector<double>(3)); // матрица жесткости
  69.         vector<vector<double>> localMatrix(3, vector<double>(3)); // локальная матрица 
  70.         vector<double> localVector(3); // локальный вектор 
  71.         this->matrixMass = matrixMass;
  72.         this->matrixStiffness = matrixStiffness;
  73.         this->localMatrix = localMatrix;
  74.         this->localVector = localVector;
  75.     }
  76. };
  77.  
  78. // краевые уловия 1 рода
  79. struct BoundaryCondition1
  80. {
  81. public:
  82.     Node node; // номер узла
  83.     int formulaNumber; // номер формулы краевого уловия 
  84.     double GetValue(double t)
  85.     {
  86.         switch (formulaNumber)
  87.         {
  88.         case 0:
  89.             return 5 * t;
  90.         case 1:
  91.             return  6 * t * node.z + 4 * t;
  92.         case 2:
  93.             return  6 * t * node.r + 4 * t;
  94.         case 3:
  95.             return 6 * t * node.z * node.z + 4 * t;
  96.         case 4:
  97.             return 6 * t * node.r * node.r + 4 * t;
  98.         case 5:
  99.             return 5 * t * t;
  100.         case 6:
  101.             return 6 * t * t * node.z + 4 * t * t;
  102.         case 7:
  103.             return 6 * t * t * node.r + 4 * t * t;
  104.         case 8:
  105.             return 6 * t * t * node.z * node.z + 4 * t * t;
  106.         case 9:
  107.             return 6 * t * t * node.r * node.r + 4 * t * t;
  108.         default:
  109.             break;
  110.         }
  111.     }
  112. };
  113.  
  114. // краевые уловия 2 рода
  115. struct BoundaryCondition2
  116. {
  117. public:
  118.     Node node[2]; // номера узлов грани
  119.     int formulaNumber; // номер формулы краевого уловия
  120.  
  121.     double GetValue(double t)
  122.     {
  123.         return b[formulaNumber](node, t);
  124.     }
  125.  
  126. private:
  127.     vector<function<double(Node node[2], double t)>> b = { [](Node node[2], double t) { return 21; } };
  128. };
  129.  
  130. struct Matrix
  131. {
  132.     vector<int> ig;
  133.     vector<int> jg;
  134.     vector<double> di;
  135.     vector<double> gg;
  136. };
  137.  
  138. void GenerationMatrixMasStifVec(Triangle& tr, double t)
  139. {
  140.     double detD = (tr.node[1].r - tr.node[0].r) * (tr.node[2].z - tr.node[0].z) - (tr.node[2].r - tr.node[0].r) * (tr.node[1].z - tr.node[0].z);
  141.  
  142.     double alphaMatrix[3][2] =
  143.     {
  144.         {(tr.node[1].z - tr.node[2].z) / detD, (tr.node[2].r - tr.node[1].r) / detD},
  145.         {(tr.node[2].z - tr.node[0].z) / detD, (tr.node[0].r - tr.node[2].r) / detD},
  146.         {(tr.node[0].z - tr.node[1].z) / detD, (tr.node[1].r - tr.node[0].r) / detD}
  147.     };
  148.  
  149.     detD = abs(detD);
  150.  
  151.     double r1[3] =
  152.     {
  153.         2 * tr.node[0].r + 2 * tr.node[1].r + tr.node[2].r,
  154.         2 * tr.node[2].r + 2 * tr.node[0].r + tr.node[1].r,
  155.         2 * tr.node[2].r + 2 * tr.node[1].r + tr.node[0].r
  156.     };
  157.  
  158.     double r2[3] =
  159.     {
  160.         3 * tr.node[0].r + tr.node[1].r + tr.node[2].r,
  161.         3 * tr.node[1].r + tr.node[0].r + tr.node[2].r,
  162.         3 * tr.node[2].r + tr.node[1].r + tr.node[0].r
  163.     };
  164.  
  165.     double temp = (tr.range.diffusionCoefficient * detD * (tr.node[0].r + tr.node[1].r + tr.node[2].r)) / 6.0;
  166.  
  167.     double temp1 = (tr.range.gamma * detD) / 120.0;
  168.     double temp2 = (tr.range.gamma * detD) / 60.0;
  169.  
  170.     double temp3 = detD / 120.0;
  171.  
  172.     double f[3]
  173.     {
  174.         tr.range.GetValueF(tr.node[0], t),
  175.         tr.range.GetValueF(tr.node[1], t),
  176.         tr.range.GetValueF(tr.node[2], t)
  177.     };
  178.  
  179.     tr.localVector[0] = (2 * r2[0] * f[0] + r1[0] * f[1] + r1[1] * f[2]) * temp3;
  180.     tr.localVector[1] = (2 * r2[1] * f[1] + r1[0] * f[0] + r1[2] * f[2]) * temp3;
  181.     tr.localVector[2] = (2 * r2[2] * f[2] + r1[1] * f[0] + r1[2] * f[1]) * temp3;
  182.  
  183.     for (int i = 0; i < 3; i++)
  184.     {
  185.         tr.matrixStiffness[i][i] = temp * (alphaMatrix[i][0] * alphaMatrix[i][0] + alphaMatrix[i][1] * alphaMatrix[i][1]);
  186.         tr.matrixMass[i][i] = temp2 * r2[i];
  187.         for (int j = 0; j < i; j++)
  188.         {
  189.             tr.matrixStiffness[i][j] = temp * (alphaMatrix[i][0] * alphaMatrix[j][0] + alphaMatrix[i][1] * alphaMatrix[j][1]);
  190.             tr.matrixStiffness[j][i] = tr.matrixStiffness[i][j];
  191.             tr.matrixMass[i][j] = temp1 * r1[+ j - 1];
  192.             tr.matrixMass[j][i] = tr.matrixMass[i][j];
  193.         }
  194.     }
  195. }
  196.  
  197. Matrix FormingPortraitMatrix()
  198. {
  199.     Matrix m;
  200.     ifstream fin;
  201.     fin.open("Portrait.txt");
  202.     int n;
  203.     fin >> n;
  204.     m.ig.resize(+ 1), m.ig.shrink_to_fit();
  205.     for (int i = 0; i < n + 1; i++)
  206.         fin >> m.ig[i];
  207.     m.jg.resize(m.ig[n]), m.jg.shrink_to_fit();
  208.     for (int i = 0; i < m.ig[n]; i++)
  209.         fin >> m.jg[i];
  210.     m.di.resize(n), m.di.shrink_to_fit();
  211.     m.gg.resize(m.ig[n]), m.gg.shrink_to_fit();
  212.     fin.close();
  213.     return m;
  214. }
  215.  
  216. Matrix Portrait(vector<Triangle>& triangles, int n)
  217. {
  218.     vector<set<int>> connection(n);
  219.  
  220.     for (auto& t : triangles)
  221.     {
  222.         for (int i = 1; i < 3; i++)
  223.             for (int j = 0; j < i; j++)
  224.             {
  225.                 auto a = t.node[i];
  226.                 auto b = t.node[j];
  227.                 if (a.globalNumber < b.globalNumber) swap(a, b);
  228.  
  229.                 connection[a.globalNumber].insert(b.globalNumber);
  230.             }
  231.     }
  232.  
  233.     Matrix globalMatrix;
  234.  
  235.     globalMatrix.ig.resize(+ 1);
  236.     int* IA = &globalMatrix.ig[0];
  237.     IA[0] = IA[1] = 0;
  238.  
  239.     for (int i = 2; i <= n; i++)
  240.     {
  241.         int col = IA[- 1];
  242.         IA[i] = col + connection[- 1].size();
  243.     }
  244.  
  245.     globalMatrix.jg.resize(IA[n]);
  246.     int* JA = &globalMatrix.jg[0];
  247.     for (int i = 1, k = 0; i < n; i++)
  248.         for (int j : connection[i])
  249.         {
  250.             JA[k] = j;
  251.             k++;
  252.         }
  253.  
  254.     globalMatrix.di.resize(n);
  255.     globalMatrix.gg.resize(IA[n]);
  256.     return globalMatrix;
  257. }
  258.  
  259. void EntryBoundaryCondition1(vector<BoundaryCondition1>& boundaryCondition1, Matrix& m, vector<double>& b, double t)
  260. {
  261.     double maxElement1 = *max_element(m.di.begin(), m.di.end());
  262.     double maxElement2 = *max_element(m.gg.begin(), m.gg.end());
  263.     double max;
  264.     (maxElement1 < maxElement2) ? max = maxElement2 * 1.0E+30 : max = maxElement1 * 1.0E+30;
  265.     for (int i = 0; i < boundaryCondition1.size(); i++)
  266.     {
  267.         m.di[boundaryCondition1[i].node.globalNumber] = max;
  268.         b[boundaryCondition1[i].node.globalNumber] = max * boundaryCondition1[i].GetValue(t);
  269.     }
  270. }
  271.  
  272. void EntryBoundaryCondition2(vector<BoundaryCondition2>& boundaryCondition2, vector<double>& b, double t)
  273. {
  274.     double h;
  275.     for (int i = 0; i < boundaryCondition2.size(); i++)
  276.     {
  277.         h = sqrt((boundaryCondition2[i].node[0].r - boundaryCondition2[i].node[1].r) * (boundaryCondition2[i].node[0].r - boundaryCondition2[i].node[1].r)
  278.             + (boundaryCondition2[i].node[0].z - boundaryCondition2[i].node[1].z) * (boundaryCondition2[i].node[0].z - boundaryCondition2[i].node[1].z));
  279.         b[boundaryCondition2[i].node[0].globalNumber] += (boundaryCondition2[i].GetValue(t) * h * (2 * boundaryCondition2[i].node[0].r + boundaryCondition2[i].node[1].r)) / 6.0;
  280.         b[boundaryCondition2[i].node[1].globalNumber] += (boundaryCondition2[i].GetValue(t) * h * (2 * boundaryCondition2[i].node[1].r + boundaryCondition2[i].node[0].r)) / 6.0;
  281.     }
  282. }
  283.  
  284. void GenerationGlobalMatrixVector(Matrix& M, vector<double>& B, vector<Triangle>& triangle)
  285. {
  286.     vector<int> L(3);
  287.     int a, b;
  288.     for (int l = 0; l < triangle.size(); l++)
  289.     {
  290.         L = { triangle[l].node[0].globalNumber, triangle[l].node[1].globalNumber, triangle[l].node[2].globalNumber };
  291.         //заносим диагональные элементы
  292.  
  293.         //начинаем цикл по строкам нижнего(и одновременно по столбцам верхнег) треугольника локальной матрицы
  294.         for (int i = 0; i < 3; i++)
  295.         {
  296.             M.di[L[i]] += triangle[l].localMatrix[i][i];
  297.             B[L[i]] += triangle[l].localVector[i];
  298.             // начинаем поиск по строке нижнего (и по столбцу верхнего) треугольника
  299.             for (int j = 0; j < i; j++)
  300.             {
  301.                 a = L[i];
  302.                 b = L[j];
  303.                 if (< b) swap(a, b);
  304.                 auto beg = M.jg.begin();
  305.                 beg += M.ig[a];
  306.                 if (M.ig[+ 1] > M.ig[a])
  307.                 {
  308.                     auto end = M.jg.begin();
  309.                     end += M.ig[+ 1] - 1;
  310.                     auto iter = lower_bound(beg, end, b);
  311.                     auto index = iter - M.jg.begin();
  312.                     M.gg[index] += triangle[l].localMatrix[i][j];
  313.                 }
  314.  
  315.             }
  316.         }
  317.  
  318.     }
  319. }
  320.  
  321. vector<vector<double>> MultMatrixNumber(vector<vector<double>>& M, double t)
  322. {
  323.     vector<vector<double>> result(M.size(), vector<double>(M.size()));
  324.     for (int i = 0; i < M.size(); i++)
  325.     {
  326.         result[i][i] = M[i][i] * t;
  327.         for (int j = 0; j < i; j++)
  328.         {
  329.             result[i][j] = M[i][j] * t;
  330.             result[j][i] = result[i][j];
  331.         }
  332.     }
  333.     return result;
  334. }
  335.  
  336. vector<vector<double>> MatrixAddition(vector<vector<double>>& A1, vector<vector<double>>& A2)
  337. {
  338.     vector<vector<double>> result(A1.size(), vector<double>(A1.size()));
  339.     for (int i = 0; i < A1.size(); i++)
  340.     {
  341.         result[i][i] = A1[i][i] + A2[i][i];
  342.         for (int j = 0; j < i; j++)
  343.         {
  344.             result[i][j] = A1[i][j] + A2[i][j];
  345.             result[j][i] = result[i][j];
  346.         }
  347.     }
  348.     return result;
  349. }
  350.  
  351. vector<double> MultMatrixVector(vector<vector<double>>& M, vector<double>& vec)
  352. {
  353.     vector<double> result(vec.size());
  354.     for (int i = 0; i < M.size(); i++)
  355.         for (int j = 0; j < vec.size(); j++)
  356.         {
  357.             result[i] += M[i][j] * vec[j];
  358.         }
  359.     return result;
  360. }
  361.  
  362. vector<double> AdditionVectors(vector<double>& vec1, vector<double>& vec2)
  363. {
  364.     vector<double> result(vec1.size());
  365.     for (int i = 0; i < vec1.size(); i++)
  366.         result[i] = vec1[i] + vec2[i];
  367.     return result;
  368. }
  369.  
  370. //вычисление нормы 
  371. double NormVector(vector<double>& x)
  372. {
  373.     double norma = 0;
  374.     for (int i = 0; i < x.size(); i++)
  375.         norma += x[i] * x[i];
  376.     norma = sqrt(norma);
  377.     return norma;
  378. }
  379.  
  380. //скалярное произведение векторов 
  381. double ScalarProduct(vector<double>& x, vector<double>& y)
  382. {
  383.     double scalarProduct = 0;
  384.     for (int i = 0; i < x.size(); i++)
  385.         scalarProduct += x[i] * y[i];
  386.     return scalarProduct;
  387. }
  388.  
  389. vector<double> MatrixVectorProduct(Matrix& A, vector<double>& x)
  390. {
  391.     vector<double> resultVector(x.size());
  392.     for (int i = 0; i < A.di.size(); i++)
  393.     {
  394.         resultVector[i] = A.di[i] * x[i];
  395.         for (int j = A.ig[i]; j < A.ig[+ 1]; j++)
  396.         {
  397.             resultVector[i] += A.gg[j] * x[A.jg[j]];
  398.             resultVector[A.jg[j]] += A.gg[j] * x[i];
  399.         }
  400.     }
  401.     return resultVector;
  402. }
  403.  
  404. vector<double> MCG(Matrix& A, vector<double>& B)
  405. {
  406.     int n = A.di.size();
  407.     vector<double> composition_Az(n);
  408.     vector<double> z = B;
  409.     vector<double> r = B;
  410.     vector<double> approximation(n);
  411.     double scalarProduct_rr, normVectorF = NormVector(B);
  412.     double relativeDiscrepancy;
  413.     double a, b, eps = 1.0E-35;
  414.     int iteration = 0, maxIteration = 100;
  415.     scalarProduct_rr = ScalarProduct(r, r);
  416.     do
  417.     {
  418.         composition_Az = MatrixVectorProduct(A, z);
  419.         a = scalarProduct_rr / ScalarProduct(composition_Az, z);
  420.         for (int i = 0; i < n; i++)
  421.             approximation[i] += a * z[i];
  422.         for (int i = 0; i < n; i++)
  423.             r[i] -= a * composition_Az[i];
  424.         b = scalarProduct_rr;
  425.         scalarProduct_rr = ScalarProduct(r, r);
  426.         if (abs(b) < 1.0e-60)
  427.             break;
  428.         b = scalarProduct_rr / b;
  429.         for (int i = 0; i < n; i++)
  430.             z[i] = z[i] * b + r[i];
  431.         relativeDiscrepancy = NormVector(r) / normVectorF;
  432.         iteration++;
  433.     } while (iteration < maxIteration);
  434.     return approximation;
  435. }
  436.  
  437. vector<double> TwoLayerScheme(Matrix& portret, vector<Triangle>& triangle, vector<BoundaryCondition1>& boundaryCondition1, vector<BoundaryCondition2>& boundaryCondition2, vector<double> qLast, vector<double>& t)
  438. {
  439.     vector<double>tmpV(3);
  440.     vector<vector<double>> tmp;
  441.     for (int i = 0; i < triangle.size(); i++)
  442.     {
  443.         tmpV[0] = qLast[triangle[i].node[0].globalNumber], tmpV[1] = qLast[triangle[i].node[1].globalNumber], tmpV[2] = qLast[triangle[i].node[2].globalNumber];
  444.         tmp = MultMatrixNumber(triangle[i].matrixMass, 1.0 / (t[1] - t[0]));
  445.         triangle[i].localMatrix = MatrixAddition(tmp, triangle[i].matrixStiffness);
  446.         tmpV = MultMatrixVector(tmp, tmpV);
  447.         triangle[i].localVector = AdditionVectors(triangle[i].localVector, tmpV);
  448.     }
  449.     Matrix M = portret;
  450.     vector<double> B(M.di.size());
  451.     GenerationGlobalMatrixVector(M, B, triangle);
  452.     EntryBoundaryCondition1(boundaryCondition1, M, B, t[1]);
  453.     EntryBoundaryCondition2(boundaryCondition2, B, t[1]);
  454.     return MCG(M, B);
  455. }
  456.  
  457. vector<double> ThreeLayerScheme(Matrix& portret, vector<Triangle>& triangle, vector<BoundaryCondition1>& boundaryCondition1, vector<BoundaryCondition2>& boundaryCondition2, vector<double> qLast, vector<double> qPreLast, vector<double>& t)
  458. {
  459.     vector<double>tmpV1(3), tmpV2(3);
  460.     vector<vector<double>> tmp;
  461.     for (int i = 0; i < triangle.size(); i++)
  462.     {
  463.         tmpV1[0] = qLast[triangle[i].node[0].globalNumber], tmpV1[1] = qLast[triangle[i].node[1].globalNumber], tmpV1[2] = qLast[triangle[i].node[2].globalNumber];
  464.         tmpV2[0] = qPreLast[triangle[i].node[0].globalNumber], tmpV2[1] = qPreLast[triangle[i].node[1].globalNumber], tmpV2[2] = qPreLast[triangle[i].node[2].globalNumber];
  465.         tmp = MultMatrixNumber(triangle[i].matrixMass, ((t[2] - t[0]) + (t[2] - t[1])) / ((t[2] - t[0]) * (t[2] - t[1])));
  466.         triangle[i].localMatrix = MatrixAddition(tmp, triangle[i].matrixStiffness);
  467.         tmp = MultMatrixNumber(triangle[i].matrixMass, ((-1) * (t[2] - t[1])) / ((t[2] - t[0]) * (t[1] - t[0])));
  468.         tmpV2 = MultMatrixVector(tmp, tmpV2);
  469.         tmp = MultMatrixNumber(triangle[i].matrixMass, (t[2] - t[0]) / ((t[1] - t[0]) * (t[2] - t[1])));
  470.         tmpV1 = MultMatrixVector(tmp, tmpV1);
  471.         triangle[i].localVector = AdditionVectors(triangle[i].localVector, tmpV1);
  472.         triangle[i].localVector = AdditionVectors(triangle[i].localVector, tmpV2);
  473.     }
  474.     Matrix M = portret;
  475.     vector<double> B(M.di.size());
  476.     GenerationGlobalMatrixVector(M, B, triangle);
  477.     EntryBoundaryCondition1(boundaryCondition1, M, B, t[2]);
  478.     EntryBoundaryCondition2(boundaryCondition2, B, t[2]);
  479.     return MCG(M, B);
  480. }
  481.  
  482. vector<Node> InputNodes()
  483. {
  484.     vector<Node> nodes;
  485.     int n;
  486.     ifstream fin;
  487.     fin.open("points.txt");
  488.     fin >> n;
  489.     nodes.resize(n), nodes.shrink_to_fit();
  490.     for (int i = 0; i < n; i++)
  491.     {
  492.         fin >> nodes[i].r >> nodes[i].z;
  493.         nodes[i].globalNumber = i;
  494.     }
  495.     fin.close();
  496.     return nodes;
  497. }
  498.  
  499. vector<Range> InputRange(int number)
  500. {
  501.     vector<Range> ranges;
  502.     ifstream fin;
  503.     fin.open("ranges.txt");
  504.     int n;
  505.     fin >> n;
  506.     ranges.resize(n), ranges.shrink_to_fit();
  507.     for (int i = 0; i < n; i++)
  508.     {
  509.         fin >> ranges[i].diffusionCoefficient >> ranges[i].gamma;
  510.         ranges[i].rangeNumber = number;
  511.     }
  512.     fin.close();
  513.     return ranges;
  514. }
  515.  
  516. vector<BoundaryCondition1> InputBoundaryCondition1(vector<Node>& nodes, int number)
  517. {
  518.     vector<BoundaryCondition1> boundaryCondition1;
  519.     ifstream fin;
  520.     fin.open("BoundaryCondition1.txt");
  521.     int n, temp1, temp2;
  522.     fin >> n;
  523.     boundaryCondition1.resize(n), boundaryCondition1.shrink_to_fit();
  524.     for (int i = 0; i < n; i++)
  525.     {
  526.         fin >> temp1;
  527.         boundaryCondition1[i].formulaNumber = number;
  528.         boundaryCondition1[i].node = nodes[temp1];
  529.     }
  530.     fin.close();
  531.     return boundaryCondition1;
  532. }
  533.  
  534. vector<BoundaryCondition2> InputBoundaryCondition2(vector<Node>& nodes)
  535. {
  536.     vector<BoundaryCondition2> boundaryCondition2;
  537.     ifstream fin;
  538.     int n, temp1, temp2;
  539.     fin.open("BoundaryCondition2.txt");
  540.     fin >> n;
  541.     boundaryCondition2.resize(n), boundaryCondition2.shrink_to_fit();
  542.     for (int i = 0; i < n; i++)
  543.     {
  544.         fin >> temp1 >> temp2 >> boundaryCondition2[i].formulaNumber;
  545.         boundaryCondition2[i].node[0] = nodes[temp1];
  546.         boundaryCondition2[i].node[1] = nodes[temp2];
  547.     }
  548.     fin.close();
  549.     return boundaryCondition2;
  550. }
  551.  
  552. double GetValueFor_q(int number, Node node, double t)
  553. {
  554.     switch (number)
  555.     {
  556.     case 0:
  557.         return 5 * t;
  558.     case 1:
  559.         return  6 * t * node.z + 4 * t;
  560.     case 2:
  561.         return  6 * t * node.r + 4 * t;
  562.     case 3:
  563.         return 6 * t * node.z * node.z + 4 *t;
  564.     case 4:
  565.         return 6 * t * node.r * node.r + 4 * t;
  566.     case 5:
  567.         return 5 * t * t;
  568.     case 6:
  569.         return 6 * t * t * node.z + 4 * t * t;
  570.     case 7:
  571.         return 6 * t * t * node.r + 4 * t * t;
  572.     case 8:
  573.         return 6 * t * t * node.z * node.z + 4 * t * t;
  574.     case 9:
  575.         return 6 * t * t * node.r * node.r + 4 * t * t;
  576.     default:
  577.         break;
  578.     }
  579. }
  580.  
  581. vector<double> Generationq0(double number, vector<Node>& nodes, double t)
  582. {
  583.     vector<double> q;
  584.     q.resize(nodes.size()), q.shrink_to_fit();
  585.     for (int i = 0; i < q.size(); i++)
  586.         q[i] = GetValueFor_q(number, nodes[i], t);
  587.     return q;
  588. }
  589.  
  590. vector<vector<double>> GenerationDecision(double number, vector<Node>& nodes, vector<double> t)
  591. {
  592.     vector<vector<double>> q;
  593.     q.resize(3), q.shrink_to_fit();
  594.     q[0].resize(nodes.size()), q.shrink_to_fit();
  595.     q[1].resize(nodes.size()), q.shrink_to_fit();
  596.     q[2].resize(nodes.size()), q.shrink_to_fit();
  597.     for (int i = 0; i < q[0].size(); i++)
  598.     {
  599.         q[0][i] = GetValueFor_q(number, nodes[i], t[0]);
  600.         q[1][i] = GetValueFor_q(number, nodes[i], t[1]);
  601.         q[2][i] = GetValueFor_q(number, nodes[i], t[2]);
  602.     }
  603.     return q;
  604. }
  605.  
  606. vector<double> InputTime()
  607. {
  608.     vector<double> t;
  609.     int n;
  610.     ifstream fin;
  611.     fin.open("time.txt");
  612.     fin >> n;
  613.     t.resize(n), t.shrink_to_fit();
  614.     for (int i = 0; i < n; i++)
  615.         fin >> t[i];
  616.     fin.close();
  617.     return t;
  618. }
  619.  
  620. void Output(vector<vector<double>>& q, vector<vector<double>>& decision, ofstream& out, int number)
  621. {
  622.     vector<string> func =
  623.     {
  624.         "5 * t",
  625.         "6 * t * z + 4 * t",
  626.         "6 * t * r + 4 * t",
  627.         "6 * t * z * z + 4 *t",
  628.         "6 * t * r * r + 4 * t",
  629.         "5 * t * t",
  630.         "6 * t * t * z + 4 * t * t",
  631.         "6 * t * t * r + 4 * t * t",
  632.         "6 * t * t * z * z + 4 * t * t",
  633.         "return 6 * t * t * r * r + 4 * t * t"
  634.     };
  635.     out << "Тест " << func[number] << endl;
  636.     out.setf(ios::scientific);
  637.     out.precision(7);
  638.     for (int j = 0; j < q[0].size(); j++)
  639.         out << q[0][j] << "\t" << q[1][j] << "\t" << q[2][j] << "\t" << decision[0][j] << "\t" << decision[1][j] << "\t" << decision[2][j] << endl;
  640.     out << endl;
  641. }
  642.  
  643. vector<Triangle> FormingArrayTriangles(vector<Range>& ranges, vector<Node>& nodes, double t)
  644. {
  645.     vector<Triangle> triangles;
  646.     ifstream fin;
  647.     fin.open("triangles.txt");
  648.     int n, numberNodes[3], numberRanges;
  649.     fin >> n;
  650.     triangles.resize(n), triangles.shrink_to_fit();
  651.     for (int i = 0; i < n; i++)
  652.     {
  653.         fin >> numberNodes[0] >> numberNodes[1] >> numberNodes[2] >> numberRanges;
  654.         Triangle tr;
  655.         tr.node[0] = nodes[numberNodes[0]], tr.node[1] = nodes[numberNodes[1]], tr.node[2] = nodes[numberNodes[2]], tr.range = ranges[numberRanges];
  656.         GenerationMatrixMasStifVec(tr, t);
  657.         triangles[i] = tr;
  658.     }
  659.     fin.close();
  660.     return triangles;
  661. }
  662.  
  663. int main()
  664. {
  665.     vector<Node> nodes;
  666.     vector<Range> ranges;
  667.     vector<Triangle> triangles;
  668.     vector<BoundaryCondition1> boundaryCondition1;
  669.     vector<BoundaryCondition2> boundaryCondition2;
  670.     vector<vector<double>> q, decision;
  671.     vector<double> t, time;
  672.     time = InputTime();
  673.     ofstream out("q.txt");
  674.     for (size_t i = 0; i < 10; i++)
  675.     {
  676.         nodes = InputNodes();
  677.         ranges = InputRange(i);
  678.         boundaryCondition1 = InputBoundaryCondition1(nodes, i);
  679.         boundaryCondition2 = InputBoundaryCondition2(nodes);
  680.         t.resize(2);
  681.         t[0] = time[0], t[1] = time[1];
  682.         triangles = FormingArrayTriangles(ranges, nodes, t[1]);
  683.         Matrix port = Portrait(triangles, nodes.size());
  684.         q.push_back(Generationq0(i, nodes, t[0]));
  685.         q.push_back(TwoLayerScheme(port, triangles, boundaryCondition1, boundaryCondition2, q[0], t));
  686.         t.resize(3);
  687.         t[0] = time[0], t[1] = time[1], t[2] = time[2];
  688.         triangles = FormingArrayTriangles(ranges, nodes, t[2]);
  689.         q.push_back(ThreeLayerScheme(port, triangles, boundaryCondition1, boundaryCondition2, q[1], q[0], t));
  690.         decision = GenerationDecision(i, nodes, t);
  691.         Output(q, decision, out, i);
  692.         q.clear();
  693.     }
  694.     out.close();
  695. }

Editor

You can edit this paste and save as new: