From 274d59e659dede795678bf62f16c148804fd4a91 Mon Sep 17 00:00:00 2001 From: ParkSuMin Date: Wed, 10 Sep 2025 21:02:45 +0300 Subject: [PATCH] =?UTF-8?q?=D0=92=D0=BE=D0=B7=D0=BC=D0=BE=D0=B6=D0=BD?= =?UTF-8?q?=D0=BE=20=D0=BA=D0=BE=D1=80=D1=80=D0=B5=D0=BA=D1=82=D0=BD=D0=BE?= =?UTF-8?q?=D0=B5=20=D1=80=D0=B5=D1=88=D0=B5=D0=BD=D0=B8=D0=B5=20=D0=BF?= =?UTF-8?q?=D1=80=D0=B8=20=D0=BF=D0=BE=D0=BC=D0=BE=D1=89=D0=B8=20=D0=BA?= =?UTF-8?q?=D1=83=D0=B1=D0=B8=D1=87=D0=B5=D1=81=D0=BA=D0=BE=D0=B3=D0=BE=20?= =?UTF-8?q?=D0=9A=D0=AD?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Solver.cpp | 119 ++++++++++++++++++++++++++++++++++++----------------- 1 file changed, 81 insertions(+), 38 deletions(-) diff --git a/Solver.cpp b/Solver.cpp index b2fe55a..95e180c 100644 --- a/Solver.cpp +++ b/Solver.cpp @@ -46,8 +46,8 @@ void Solver::Execute_Linear() { std::cout << "Ansamb load vector:\n" << global_load << std::endl; // ГУ 1-го рода - double u_left = 5.0; // u(1) = 5 - double u_right = 15.0; // u(6) = 15 + double u_left = -5.0; // u(1) = 5 + double u_right = -10.0; // u(6) = 15 // 1. Учёт условия u(1)=5 for (int i = 1; i < N + 1; ++i) { @@ -85,56 +85,99 @@ void Solver::Execute_Linear() { } void Solver::Execute_Cubic() { + // Локальная матрица жесткости (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 / 80; - local(1, 0) = 189/(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) = -20 * 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, 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, 3) = 189 * A / (40 * dx) + 57 * B / 80; - local(3, 0) = 13*A/(40 * dx); + 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, 2) = 189 * A / (40 * dx) - 57 * B / 80; + local(3, 3) = -37 * A / (10 * dx) + B / 2; + // Локальный вектор нагрузки (4x1) VectorXd local_load(4); - local_load(0) = local_load(3) = -C * dx / 8; - local_load(1) = local_load(2) = -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; - // Ансамблированные матрица и вектор - MatrixXd ansamb = MatrixXd::Zero(N + 1, N + 1); - VectorXd global_load = VectorXd::Zero(N + 1); + // Размер глобальной системы: для кубических элементов (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_i = elem; - // int node_j = elem + 1; + // АНСАМБЛИРОВАНИЕ + for (int elem = 0; elem < N; ++elem) { + int node_start = 2 * elem; // Начальный индекс для текущего элемента - // // Матрица жесткости - // 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); + // Добавляем локальную матрицу + for (int i = 0; i < 4; ++i) { + for (int j = 0; j < 4; ++j) { + ansamb(node_start + i, node_start + j) += local(i, j); + } + } - // // Вектор нагрузки - // global_load(node_i) += local_load(0); - // global_load(node_j) += local_load(1); - //} + // Добавляем локальный вектор нагрузки + for (int i = 0; i < 4; ++i) { + global_load(node_start + i) += local_load(i); + } + } - //std::cout << "Before:" << std::endl; - //std::cout << "Ansamb matrix:\n" << ansamb << std::endl; - //std::cout << "Ansamb load vector:\n" << global_load << std::endl; + std::cout << "Before boundary conditions:" << std::endl; + 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 + + // Левая граница + 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; + global_load(ndof - 1) = u_right; + + std::cout << "\nAfter boundary conditions:" << std::endl; + 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) + std::ofstream file("matrix_cubic.txt"); + for (int i = 0; i < ndof; i += 2) { + file << solution(i) << ' '; + } + file << std::endl; } \ No newline at end of file