diff --git a/MemMAPR-MKE.cpp b/MemMAPR-MKE.cpp index 40a126d..e20a88d 100644 --- a/MemMAPR-MKE.cpp +++ b/MemMAPR-MKE.cpp @@ -2,8 +2,10 @@ #include "Solver.h" int main() { - Solver slv(1., 0., 1., 20, 0, 1); - slv.Execute_Linear(); - + Solver slv(3., 2., 5., 20, 1, 6); + std::cout << "Linear element:" << std::endl; + slv.Execute_Linear(-5, -10); + std::cout << "\nCubic element:" << std::endl; + slv.Execute_Cubic(-5, -10); return 0; } \ No newline at end of file diff --git a/Solver.cpp b/Solver.cpp index 95e180c..e96d9c8 100644 --- a/Solver.cpp +++ b/Solver.cpp @@ -1,4 +1,4 @@ -#include "Header.h" +#include "Header.h" #include using namespace Eigen; @@ -9,34 +9,34 @@ Solver::Solver(double _A, double _B, double _C, int _N, int _l, int _u) { dx = L / N; } -void Solver::Execute_Linear() { +void Solver::Execute_Linear(double val1, double val2) { MatrixXd local = MatrixXd::Zero(2, 2); - local(0, 0) = A / dx - B / 2; - local(0, 1) = -A / dx + B / 2; - local(1, 0) = -A / dx - B / 2; - local(1, 1) = A / dx + B / 2; + local(0, 0) = A / dx - B / 2.; + local(0, 1) = -A / dx + B / 2.; + local(1, 0) = -A / dx - B / 2.; + local(1, 1) = A / dx + B / 2.; std::cout << "Local matrix:\n" << local << std::endl; VectorXd local_load(2); local_load(0) = C * dx / 2; local_load(1) = C * dx / 2; - // + // Ансаблированные матрицы и вектор MatrixXd ansamb = MatrixXd::Zero(N + 1, N + 1); VectorXd global_load = VectorXd::Zero(N + 1); - // + // Ансамблирование для каждого конечного элемента for (int elem = 0; elem < N; ++elem) { int node_i = elem; int node_j = elem + 1; - // + // Матрица жесткости ansamb(node_i, node_i) += local(0, 0); ansamb(node_i, node_j) += local(0, 1); ansamb(node_j, node_i) += local(1, 0); ansamb(node_j, node_j) += local(1, 1); - // + // Вектор нагрузки global_load(node_i) += local_load(0); global_load(node_j) += local_load(1); } @@ -45,25 +45,17 @@ void Solver::Execute_Linear() { std::cout << "Ansamb matrix:\n" << ansamb << std::endl; std::cout << "Ansamb load vector:\n" << global_load << std::endl; - // 1- - double u_left = -5.0; // u(1) = 5 - double u_right = -10.0; // u(6) = 15 + // ГУ 1-го рода + double u_left = val1; // u(1) = 5 + double u_right = val2; // u(6) = 15 - // 1. u(1)=5 - for (int i = 1; i < N + 1; ++i) { - global_load(i) -= ansamb(i, 0) * u_left; - } - // , = 1 + // Обнуляем первую строку и столбец матрицы, устанавливаем диагональ = 1 ansamb.row(0).setZero(); ansamb.col(0).setZero(); ansamb(0, 0) = 1.0; global_load(0) = u_left; - // 2. u(6)=15 - for (int i = 0; i < N; ++i) { - global_load(i) -= ansamb(i, N) * u_right; - } - // , = 1 + // Обнуляем последнюю строку и столбец матрицы, устанавливаем диагональ = 1 ansamb.row(N).setZero(); ansamb.col(N).setZero(); ansamb(N, N) = 1.0; @@ -84,56 +76,56 @@ void Solver::Execute_Linear() { file << std::endl; } -void Solver::Execute_Cubic() { - // (4x4) +void Solver::Execute_Cubic(double val1, double val2) { + // Локальная матрица жесткости (4x4) MatrixXd local = MatrixXd::Zero(4, 4); - local(0, 0) = -37 * A / (10 * dx) - B / 2; - local(0, 1) = 189 * A / (40 * dx) + 57 * B / 80; - local(0, 2) = -27 * A / (20 * dx) - 3 * B / 10; - local(0, 3) = 13 * A / (40 * dx) + 7 * B / 80; + local(0, 0) = -37. * A / (10 * dx) - B / 2.; + local(0, 1) = 189. * A / (40 * dx) + 57. * B / 80.; + local(0, 2) = -27. * A / (20 * dx) - 3. * B / 10.; + local(0, 3) = 13. * A / (40. * dx) + 7. * B / 10.; - local(1, 0) = 189 * A / (40 * dx) - 57 * B / 80; - local(1, 1) = -54 * A / (5 * dx); - local(1, 2) = 297 * A / (40 * dx) + 81 * B / 80; - local(1, 3) = -20 * A / (20 * dx) - 3 * B / 10; + local(1, 0) = 189. * A / (40 * dx) - 57. * B / 80.; + local(1, 1) = -54. * A / (5 * dx); + local(1, 2) = 297. * A / (40. * dx) + 81. * B / 80.; + local(1, 3) = -27 * A / (20. * dx) - 3. * B / 10.; - local(2, 0) = -27 * A / (20 * dx) + 3 * B / 10; - local(2, 1) = 297 * A / (40 * dx) - 81 * B / 80; - local(2, 2) = -54 * A / (5 * dx); - local(2, 3) = 189 * A / (40 * dx) + 57 * B / 80; + local(2, 0) = -27. * A / (20. * dx) + 3. * B / 10.; + local(2, 1) = 297. * A / (40. * dx) - 81. * B / 80.; + local(2, 2) = -54. * A / (5. * dx); + local(2, 3) = 189. * A / (40. * dx) + 57. * B / 80.; - local(3, 0) = 13 * A / (40 * dx); - local(3, 1) = -27 * A / (20 * dx) + 3 * B / 10; - local(3, 2) = 189 * A / (40 * dx) - 57 * B / 80; - local(3, 3) = -37 * A / (10 * dx) + B / 2; + local(3, 0) = 13. * A / (40. * dx) - 7. * B / 80.; + local(3, 1) = -27. * A / (20. * dx) + 3. * B / 10.; + local(3, 2) = 189. * A / (40. * dx) - 57. * B / 80.; + local(3, 3) = -37. * A / (10. * dx) + B / 2.; - // (4x1) + // Локальный вектор нагрузки (4x1) VectorXd local_load(4); - local_load(0) = -C * dx / 8; - local_load(1) = -3 * C * dx / 8; - local_load(2) = -3 * C * dx / 8; - local_load(3) = -C * dx / 8; + local_load(0) = -C * dx / 8.; + local_load(1) = -3. * C * dx / 8.; + local_load(2) = -3. * C * dx / 8.; + local_load(3) = -C * dx / 8.; std::cout << "Local matrix:\n" << local << std::endl; std::cout << "Local load vector:\n" << local_load << std::endl; - // : (4 , 2 ) + // Размер глобальной системы: для кубических элементов (4 узла на элемент, перекрытие по 2 узла) int ndof = 2 * N + 2; MatrixXd ansamb = MatrixXd::Zero(ndof, ndof); VectorXd global_load = VectorXd::Zero(ndof); - // + // АНСАМБЛИРОВАНИЕ for (int elem = 0; elem < N; ++elem) { - int node_start = 2 * elem; // + int node_start = 2 * elem; // Начальный индекс для текущего элемента - // + // Добавляем локальную матрицу for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) { ansamb(node_start + i, node_start + j) += local(i, j); } } - // + // Добавляем локальный вектор нагрузки for (int i = 0; i < 4; ++i) { global_load(node_start + i) += local_load(i); } @@ -143,23 +135,15 @@ void Solver::Execute_Cubic() { std::cout << "Ansamb matrix:\n" << ansamb << std::endl; std::cout << "Ansamb load vector:\n" << global_load << std::endl; - // - double u_left = -5.0; // u(1) = 5 - double u_right = -10.0; // u(6) = 15 + // ГРАНИЧНЫЕ УСЛОВИЯ + double u_left = val1; // u(1) = 5 + double u_right = val2; // u(6) = 15 - // - for (int i = 1; i < ndof; ++i) { - global_load(i) -= ansamb(i, 0) * u_left; - } ansamb.row(0).setZero(); ansamb.col(0).setZero(); ansamb(0, 0) = 1.0; global_load(0) = u_left; - // - for (int i = 0; i < ndof - 1; ++i) { - global_load(i) -= ansamb(i, ndof - 1) * u_right; - } ansamb.row(ndof - 1).setZero(); ansamb.col(ndof - 1).setZero(); ansamb(ndof - 1, ndof - 1) = 1.0; @@ -169,12 +153,12 @@ void Solver::Execute_Cubic() { std::cout << "Modified matrix:\n" << ansamb << std::endl; std::cout << "Modified load vector:\n" << global_load << std::endl; - // + // Решение системы VectorXd solution = ansamb.fullPivLu().solve(global_load); std::cout << "\nSolution:" << std::endl; std::cout << solution << std::endl; - // ( , 2) + // Сохранение результатов (берем только значения функции в узлах, шаг 2) std::ofstream file("matrix_cubic.txt"); for (int i = 0; i < ndof; i += 2) { file << solution(i) << ' '; diff --git a/Solver.h b/Solver.h index 6eb0c8d..94effc1 100644 --- a/Solver.h +++ b/Solver.h @@ -5,6 +5,6 @@ private: int N, upper, lower; public: Solver(double _A, double _B, double _C, int _N, int _l, int _u); - void Execute_Linear(); - void Execute_Cubic(); + void Execute_Linear(double value_1, double value_2); + void Execute_Cubic(double value_1, double value_2); }; \ No newline at end of file