ParkSuMin/MemMAPR-MKE
1
0
Fork 1
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Complete cubic solver

Browse Source 
Изменены коэффициенты в локальных матрицах жёсткости. Значения int приведены к double формату
This commit is contained in:
ParkSuMin
2025-09-12 23:23:26 +03:00
parent 274d59e659
commit e53cef771a
3 changed files with 55 additions and 69 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
MemMAPR-MKE.cpp
View File

@@ -2,8 +2,10 @@
#include "Solver.h" #include "Solver.h"
int main() { int main() {
Solver slv(1., 0., 1., 20, 0, 1); Solver slv(3., 2., 5., 20, 1, 6);
slv.Execute_Linear(); std::cout << "Linear element:" << std::endl;
slv.Execute_Linear(-5, -10);
std::cout << "\nCubic element:" << std::endl;
slv.Execute_Cubic(-5, -10);
return 0; return 0;
} }

112
Solver.cpp
View File

@@ -1,4 +1,4 @@
#include "Header.h" #include "Header.h"
#include <Eigen/Dense> #include <Eigen/Dense>
using namespace Eigen; using namespace Eigen;
@@ -9,34 +9,34 @@ Solver::Solver(double _A, double _B, double _C, int _N, int _l, int _u) {
dx = L / N; dx = L / N;
} }
void Solver::Execute_Linear() { void Solver::Execute_Linear(double val1, double val2) {
MatrixXd local = MatrixXd::Zero(2, 2); MatrixXd local = MatrixXd::Zero(2, 2);
local(0, 0) = A / dx - B / 2; local(0, 0) = A / dx - B / 2.;
local(0, 1) = -A / dx + B / 2; local(0, 1) = -A / dx + B / 2.;
local(1, 0) = -A / dx - B / 2; local(1, 0) = -A / dx - B / 2.;
local(1, 1) = A / dx + B / 2; local(1, 1) = A / dx + B / 2.;
std::cout << "Local matrix:\n" << local << std::endl; std::cout << "Local matrix:\n" << local << std::endl;
VectorXd local_load(2); VectorXd local_load(2);
local_load(0) = C * dx / 2; local_load(0) = C * dx / 2;
local_load(1) = C * dx / 2; local_load(1) = C * dx / 2;
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> // Ансаблированные матрицы и вектор
MatrixXd ansamb = MatrixXd::Zero(N + 1, N + 1); MatrixXd ansamb = MatrixXd::Zero(N + 1, N + 1);
VectorXd global_load = VectorXd::Zero(N + 1); VectorXd global_load = VectorXd::Zero(N + 1);
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> // Ансамблирование для каждого конечного элемента
for (int elem = 0; elem < N; ++elem) { for (int elem = 0; elem < N; ++elem) {
int node_i = elem; int node_i = elem;
int node_j = elem + 1; int node_j = elem + 1;
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> // Матрица жесткости
ansamb(node_i, node_i) += local(0, 0); ansamb(node_i, node_i) += local(0, 0);
ansamb(node_i, node_j) += local(0, 1); ansamb(node_i, node_j) += local(0, 1);
ansamb(node_j, node_i) += local(1, 0); ansamb(node_j, node_i) += local(1, 0);
ansamb(node_j, node_j) += local(1, 1); ansamb(node_j, node_j) += local(1, 1);
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> // Вектор нагрузки
global_load(node_i) += local_load(0); global_load(node_i) += local_load(0);
global_load(node_j) += local_load(1); 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 matrix:\n" << ansamb << std::endl;
std::cout << "Ansamb load vector:\n" << global_load << std::endl; std::cout << "Ansamb load vector:\n" << global_load << std::endl;
// <EFBFBD><EFBFBD> 1-<EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD> // ГУ 1-го рода
double u_left = -5.0; // u(1) = 5 double u_left = val1; // u(1) = 5
double u_right = -10.0; // u(6) = 15 double u_right = val2; // u(6) = 15
// 1. <20><><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> u(1)=5 // Обнуляем первую строку и столбец матрицы, устанавливаем диагональ = 1
for (int i = 1; i < N + 1; ++i) {
global_load(i) -= ansamb(i, 0) * u_left;
}
// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>, <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> = 1
ansamb.row(0).setZero(); ansamb.row(0).setZero();
ansamb.col(0).setZero(); ansamb.col(0).setZero();
ansamb(0, 0) = 1.0; ansamb(0, 0) = 1.0;
global_load(0) = u_left; global_load(0) = u_left;
// 2. <20><><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> u(6)=15 // Обнуляем последнюю строку и столбец матрицы, устанавливаем диагональ = 1
for (int i = 0; i < N; ++i) {
global_load(i) -= ansamb(i, N) * u_right;
}
// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>, <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> = 1
ansamb.row(N).setZero(); ansamb.row(N).setZero();
ansamb.col(N).setZero(); ansamb.col(N).setZero();
ansamb(N, N) = 1.0; ansamb(N, N) = 1.0;
@@ -84,56 +76,56 @@ void Solver::Execute_Linear() {
file << std::endl; file << std::endl;
} }
void Solver::Execute_Cubic() { void Solver::Execute_Cubic(double val1, double val2) {
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> (4x4) // Локальная матрица жесткости (4x4)
MatrixXd local = MatrixXd::Zero(4, 4); MatrixXd local = MatrixXd::Zero(4, 4);
local(0, 0) = -37 * A / (10 * dx) - B / 2; local(0, 0) = -37. * A / (10 * dx) - B / 2.;
local(0, 1) = 189 * A / (40 * dx) + 57 * B / 80; local(0, 1) = 189. * A / (40 * dx) + 57. * B / 80.;
local(0, 2) = -27 * A / (20 * dx) - 3 * B / 10; local(0, 2) = -27. * A / (20 * dx) - 3. * B / 10.;
local(0, 3) = 13 * A / (40 * dx) + 7 * B / 80; local(0, 3) = 13. * A / (40. * dx) + 7. * B / 10.;
local(1, 0) = 189 * A / (40 * dx) - 57 * B / 80; local(1, 0) = 189. * A / (40 * dx) - 57. * B / 80.;
local(1, 1) = -54 * A / (5 * dx); local(1, 1) = -54. * A / (5 * dx);
local(1, 2) = 297 * A / (40 * dx) + 81 * B / 80; local(1, 2) = 297. * A / (40. * dx) + 81. * B / 80.;
local(1, 3) = -20 * A / (20 * dx) - 3 * B / 10; local(1, 3) = -27 * A / (20. * dx) - 3. * B / 10.;
local(2, 0) = -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, 1) = 297. * A / (40. * dx) - 81. * B / 80.;
local(2, 2) = -54 * A / (5 * dx); 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) - 7. * B / 80.;
local(3, 1) = -27 * A / (20 * dx) + 3 * B / 10; local(3, 1) = -27. * A / (20. * dx) + 3. * B / 10.;
local(3, 2) = 189 * A / (40 * dx) - 57 * B / 80; local(3, 2) = 189. * A / (40. * dx) - 57. * B / 80.;
local(3, 3) = -37 * A / (10 * dx) + B / 2; local(3, 3) = -37. * A / (10. * dx) + B / 2.;
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> (4x1) // Локальный вектор нагрузки (4x1)
VectorXd local_load(4); VectorXd local_load(4);
local_load(0) = -C * dx / 8; local_load(0) = -C * dx / 8.;
local_load(1) = -3 * C * dx / 8; local_load(1) = -3. * C * dx / 8.;
local_load(2) = -3 * C * dx / 8; local_load(2) = -3. * C * dx / 8.;
local_load(3) = -C * dx / 8; local_load(3) = -C * dx / 8.;
std::cout << "Local matrix:\n" << local << std::endl; std::cout << "Local matrix:\n" << local << std::endl;
std::cout << "Local load vector:\n" << local_load << std::endl; std::cout << "Local load vector:\n" << local_load << std::endl;
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>: <20><><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> (4 <20><><EFBFBD><EFBFBD> <20><> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>, <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><> 2 <20><><EFBFBD><EFBFBD>) // Размер глобальной системы: для кубических элементов (4 узла на элемент, перекрытие по 2 узла)
int ndof = 2 * N + 2; int ndof = 2 * N + 2;
MatrixXd ansamb = MatrixXd::Zero(ndof, ndof); MatrixXd ansamb = MatrixXd::Zero(ndof, ndof);
VectorXd global_load = VectorXd::Zero(ndof); VectorXd global_load = VectorXd::Zero(ndof);
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> // АНСАМБЛИРОВАНИЕ
for (int elem = 0; elem < N; ++elem) { for (int elem = 0; elem < N; ++elem) {
int node_start = 2 * elem; // <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> int node_start = 2 * elem; // Начальный индекс для текущего элемента
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> // Добавляем локальную матрицу
for (int i = 0; i < 4; ++i) { for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) { for (int j = 0; j < 4; ++j) {
ansamb(node_start + i, node_start + j) += local(i, j); ansamb(node_start + i, node_start + j) += local(i, j);
} }
} }
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> // Добавляем локальный вектор нагрузки
for (int i = 0; i < 4; ++i) { for (int i = 0; i < 4; ++i) {
global_load(node_start + i) += local_load(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 matrix:\n" << ansamb << std::endl;
std::cout << "Ansamb load vector:\n" << global_load << std::endl; std::cout << "Ansamb load vector:\n" << global_load << std::endl;
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> // ГРАНИЧНЫЕ УСЛОВИЯ
double u_left = -5.0; // u(1) = 5 double u_left = val1; // u(1) = 5
double u_right = -10.0; // u(6) = 15 double u_right = val2; // u(6) = 15
// <20><><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
for (int i = 1; i < ndof; ++i) {
global_load(i) -= ansamb(i, 0) * u_left;
}
ansamb.row(0).setZero(); ansamb.row(0).setZero();
ansamb.col(0).setZero(); ansamb.col(0).setZero();
ansamb(0, 0) = 1.0; ansamb(0, 0) = 1.0;
global_load(0) = u_left; global_load(0) = u_left;
// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
for (int i = 0; i < ndof - 1; ++i) {
global_load(i) -= ansamb(i, ndof - 1) * u_right;
}
ansamb.row(ndof - 1).setZero(); ansamb.row(ndof - 1).setZero();
ansamb.col(ndof - 1).setZero(); ansamb.col(ndof - 1).setZero();
ansamb(ndof - 1, ndof - 1) = 1.0; 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 matrix:\n" << ansamb << std::endl;
std::cout << "Modified load vector:\n" << global_load << std::endl; std::cout << "Modified load vector:\n" << global_load << std::endl;
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> // Решение системы
VectorXd solution = ansamb.fullPivLu().solve(global_load); VectorXd solution = ansamb.fullPivLu().solve(global_load);
std::cout << "\nSolution:" << std::endl; std::cout << "\nSolution:" << std::endl;
std::cout << solution << std::endl; std::cout << solution << std::endl;
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> (<28><><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20> <20><><EFBFBD><EFBFBD><EFBFBD>, <20><><EFBFBD> 2) // Сохранение результатов (берем только значения функции в узлах, шаг 2)
std::ofstream file("matrix_cubic.txt"); std::ofstream file("matrix_cubic.txt");
for (int i = 0; i < ndof; i += 2) { for (int i = 0; i < ndof; i += 2) {
file << solution(i) << ' '; file << solution(i) << ' ';

4
Solver.h
View File

@@ -5,6 +5,6 @@ private:
int N, upper, lower; int N, upper, lower;
public: public:
Solver(double _A, double _B, double _C, int _N, int _l, int _u); Solver(double _A, double _B, double _C, int _N, int _l, int _u);
void Execute_Linear(); void Execute_Linear(double value_1, double value_2);
void Execute_Cubic(); void Execute_Cubic(double value_1, double value_2);
}; };