From 85cd9c0dd066bb09063fae39e02d8dbfe2e17b45 Mon Sep 17 00:00:00 2001
From: ParkSuMin <kamalovap@paksurion.ru>
Date: Wed, 10 Sep 2025 18:51:26 +0300
Subject: [PATCH] Write cubic solver

---
 Solver.cpp | 61 +++++++++++++++++++++++++++++++++++++++++++++---------
 1 file changed, 51 insertions(+), 10 deletions(-)

diff --git a/Solver.cpp b/Solver.cpp
index ce0597e..b2fe55a 100644
--- a/Solver.cpp
+++ b/Solver.cpp
@@ -21,7 +21,7 @@ void Solver::Execute_Linear() {
     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);
 
@@ -85,15 +85,56 @@ void Solver::Execute_Linear() {
 }
 
 void Solver::Execute_Cubic() {
-    MatrixXd local = MatrixXd::Zero(2, 2);
-    local(0, 0) = A * 37/(10 * dx) + B / 2;
-    local(0, 1) = -A / dx + B / 2;
-    local(0, 1) = A;
-    local(1, 0) = -A / dx - B / 2;
-    local(1, 1) = A / dx + B / 2;
+    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(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(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;
+
+
+    VectorXd local_load(4);
+    local_load(0) = local_load(3) = -C * dx / 8;
+    local_load(1) = local_load(2) = -3 * C * dx / 8;
+
     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);
+    //}
+
+    //std::cout << "Before:" << std::endl;
+    //std::cout << "Ansamb matrix:\n" << ansamb << std::endl;
+    //std::cout << "Ansamb load vector:\n" << global_load << std::endl;
+
 }
\ No newline at end of file