From 4ab15cb13b252834b1b1ec756b5639fad7dd9c93 Mon Sep 17 00:00:00 2001 From: ParkSuMin Date: Sun, 21 Sep 2025 22:28:44 +0300 Subject: [PATCH] Ready linear element eval --- Solver.cpp | 55 +++++++++++++++++++++++++++--------------------------- 1 file changed, 28 insertions(+), 27 deletions(-) diff --git a/Solver.cpp b/Solver.cpp index 63845f1..d82257d 100644 --- a/Solver.cpp +++ b/Solver.cpp @@ -10,67 +10,68 @@ Solver::Solver(double _A, double _B, double _C, double _D, int _N, int _l, int _ } 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.; - std::cout << "Local matrix:\n" << local << std::endl; - VectorXd local_load(2); - local_load(0) = C * dx / 2; - local_load(1) = C * dx / 2; - // Ансаблированные матрицы и вектор + local(0, 0) = -A / dx - B / 2. + C * dx / 3.; + local(0, 1) = A / dx + B / 2. + C * dx / 6.; + local(1, 0) = A / dx - B / 2. + C * dx / 6.; + local(1, 1) = -A / dx + B / 2. + C * dx / 3.; + + local_load(0) = -D * dx / 2.; + local_load(1) = -D * 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); } +#if DEBUG std::cout << "Before:" << std::endl; std::cout << "Ansamb matrix:\n" << ansamb << std::endl; std::cout << "Ansamb load vector:\n" << global_load << std::endl; +#endif - // ГУ 1-го рода - double u_left = val1; // u(1) = 5 - double u_right = val2; // u(6) = 15 + double u_right = val2; - // Обнуляем первую строку и столбец матрицы, устанавливаем диагональ = 1 + // Clear first and last rows ansamb.row(0).setZero(); - ansamb.col(0).setZero(); - ansamb(0, 0) = 1.0; - global_load(0) = u_left; - - // Обнуляем последнюю строку и столбец матрицы, устанавливаем диагональ = 1 ansamb.row(N).setZero(); - ansamb.col(N).setZero(); - ansamb(N, N) = 1.0; + + // u'(0) = u(0) + ansamb(0, 0) = dx + 1; + ansamb(0, 1) = -1; + global_load(0) = 0; + + // u(10) = u_right = 5 + ansamb(N, N) = 1; global_load(N) = u_right; +#if DEBUG std::cout << "\nAfter:" << std::endl; std::cout << "Modified matrix:\n" << ansamb << std::endl; std::cout << "Modified load vector:\n" << global_load << std::endl; +#endif VectorXd solution = ansamb.fullPivLu().solve(global_load); std::cout << "\nSolution:" << std::endl; std::cout << solution << std::endl; - std::ofstream file("matrix_2.txt"); - for (int i = 0; i < N + 1; i++) { + std::ofstream file("matrix_linear.txt"); + for (int i = 0; i < N; i++) { file << solution(i) << ' '; } file << std::endl; @@ -119,7 +120,7 @@ void Solver::Execute_Cubic(double val1, double val2) { Amat(i + 3, i + 1) += B * 3.0 / 10.0; Amat(i + 3, i + 2) += B * (-57.0) / 80.0; Amat(i + 3, i + 3) += B * 1.0 / 2.0; - } + } // Assembdxe vector for (int i = 0; i < mat_dim - 3; i += 3) { @@ -143,7 +144,7 @@ void Solver::Execute_Cubic(double val1, double val2) { VectorXd solution = Amat.colPivHouseholderQr().solve(b); std::cout << "\nSolution:" << std::endl; std::cout << solution << std::endl; - + std::ofstream file("matrix_cubic.txt"); for (int i = 0; i < solution.size(); i++) { file << solution(i) << ' ';