#define _CRT_SECURE_NO_WARNINGS
#include <clocale>
#include <fstream>
#include <iostream>

// для красивого вида формул ЛЛт
#define MatrixL global_matrix_A
#define VectX vector_B
#define VectY vector_B

using namespace std;
const double lamda = 1; // лямда
const double gamma = 1; // гамма
double global_matrix_A;
double vector_B;
double grid; // указатель на массив узлов сетки
/////Для ллтГауса НАЧ/////
double ini_vect_y(int Razm, double vector_B, double global_matrix_A)
{
for (int i = 0; i < Razm; i++)
{
double sum = 0;
for (int j = 0; j < i; j++)
{
sum += MatrixL[j] VectY[j];
}
VectY = (vector_B - sum) / global_matrix_A;
}

return vector_B;
}

double ini_vect_x(int Razm, double vector_B, double global_matrix_A)
{
for (int i = Razm - 1; i >= 0; i--)
{
double sum = 0;
for (int k = i + 1; k < Razm; k++)
{
sum += MatrixL[k] VectX[k];
}
VectX = (VectY - sum) / MatrixL;
}
return vector_B;
}

double chol_decomposition(int Razm, double global_matrix_A)
{
for (int i = 0; i < Razm; i++)
{
for (int j = 0; j < i; j++)
{
double sum = 0;
for (int k = 0; k < j; k++)
{
sum += MatrixL[k] MatrixL[j][k];
}
MatrixL[j] = (global_matrix_A[j] - sum) / MatrixL[j][j];
}

double sum = 0;
for (int jj = 0; jj < i; jj++)
{
sum += MatrixL[jj] MatrixL[jj];
}
MatrixL = sqrt(global_matrix_A - sum);
}
return global_matrix_A;
}

bool sim_matrix(double global_matrix_A, int Razm)
{
int simKoof = 0; // коофицент симметричности
for (int i = 0; i < Razm; i++)
{
for (int j = 0; j < i; j++)
{
if (global_matrix_A[j] == global_matrix_A[j])
{
simKoof++;
}
}
}
if (simKoof == (Razm Razm - Razm) / 2)
{
return true;
}
else
{
cout << "Матрица не симметрична" << endl;
//system("pause");
return false;
}
}

void out_res(double vector_B, int Razm)
{
ofstream fRes;
fRes.open("Result.txt");
fRes.precision(15);
for (int i = 0; i < Razm; i++)
{
fRes << VectX << endl;
}
fRes.close();

//FILE f = fopen("Result.txt", "w");
//for (int i = 0; i < Razm; i++) fprintf(f, "%.15lf \n", VectX);
//fclose(f);
}
/////Для ллтГауса КОН/////

double f_function(double x)
{
return 2 + x;
}

double read_grid(double grid, int& quantity_of_finite_elemets, int& quantity_of_units)
{
fstream file_grid; // файл с сеткой
file_grid.open("grid.txt"); // открыть файл

if (file_grid.is_open()) // если файл открыт
{
file_grid >> quantity_of_finite_elemets; //считать количество кон элемов

quantity_of_units = quantity_of_finite_elemets + 1; //вычислить количество узлов в сетке

grid = new double[quantity_of_units]; // выделить память под сетку

for (int i = 0; i < quantity_of_units; i++) // считать узлы в сетку
{
file_grid >> grid;
}

file_grid.close();
}
else // иначе сообщение об ошибке
{
cout << "Файл не удалось открыть";
}
return grid;
}

void read_boundary_cons(int& type_of_left_boundary_con, int& type_of_right_boundary_con)
{
fstream file_kr_usl; // файл с краевыми условиями
file_kr_usl.open("types_of_boundary_cons.txt"); // открыть файл
if (file_kr_usl.is_open()) // если файл открыт
{
file_kr_usl >> type_of_left_boundary_con; //считать количество кон элемов
file_kr_usl >> type_of_right_boundary_con;

}
else // иначе сообщение об ошибке
{
cout << "Файл с краевыми условиями не удалось открыть";
}

}

double create_local_stiffnes_matrix(double lamda, double h)
{
double local_stiffnes_matrix = new double[2];
for (int i = 0; i < 2; i++)
{
local_stiffnes_matrix = new double[2];
}

// наверное стоит сделать так чтобы функция фозвращала массив с локальной матрицей
// в принципе можно передавать указаель по ссылке или делать через ретерн
// вычисление локальной матрицы жесткости для i кон элема
for (int i = 0; i < 2; i++)
{
for (int j = 0; j < 2; j++)
{
if (i == j)
{
local_stiffnes_matrix[j] = lamda / h;
}
else
{
local_stiffnes_matrix[j] = -lamda / h;
}
}
}
return local_stiffnes_matrix;
}

double create_local_mass_matrix(double gamma, double h)
{
double local_mass_matrix = new double[2];
for (int i = 0; i < 2; i++)
{
local_mass_matrix = new double[2];
}

for (int i = 0; i < 2; i++)
{
for (int j = 0; j < 2; j++)
{
if (i == j)
{
local_mass_matrix[j] = 2 (gamma h) / 6;
}
else
{
local_mass_matrix[j] = (gamma h) / 6;
}
}
}
return local_mass_matrix;
}

void set_left_boundary_con_1()
{
global_matrix_A[0][0] = 1E+050;
vector_B[0] = vector_B[0] 1E+050;

}

void set_right_boundary_con_1(int quantity_of_units)
{
global_matrix_A[quantity_of_units][quantity_of_units] = 1E+60;
vector_B[0] = vector_B[quantity_of_units] 1E+60;
}

void output() {

}

void global_compilator(int quantity_of_units, int quantity_of_finite_elements, double grid)
{
// жертвуя памятью приобретем простоту(используем плотный формат)
//выделение памяти под глобальную матрицу
double local_stiffnes_matrix = new double[2];
for (int i = 0; i < 2; i++)
{
local_stiffnes_matrix = new double[2];
}

double local_mass_matrix = new double[2];
for (int i = 0; i < 2; i++)
{
local_mass_matrix = new double[2];
}

for (int i = 0; i < quantity_of_finite_elements; i++)
{
double h = grid[i + 1] - grid;
local_stiffnes_matrix = create_local_stiffnes_matrix(lamda, h);
local_mass_matrix = create_local_mass_matrix(gamma, h);
for (int c = i; c < i + 2; c++)
{
for (int j = i; j < i + 2; j++)
{
global_matrix_A[c][j] += local_stiffnes_matrix[c - i][j - i];
global_matrix_A[c][j] += local_mass_matrix[c - i][j - i];

}
}

vector_B += h / 6 (2 f_function(grid) + f_function(grid[i + 1]));
vector_B[i + 1] += h / 6 (f_function(grid) + 2 f_function(grid[i + 1]));

}

}
void main()
{
setlocale(LC_ALL, "Russian");
// рассчетная область
// краевые условия 1 и 2 род на левом конце на правом конце
// количество конечных элементов. Зачем это нужно? Возьму длину расчетной области поделю на кол-во к элемов
// таким образом получу шаг. Этот шаг будет определять длину каждого начального элемента
// как формировать сетку?
/
Расчетная область:
Она записана в файле input/uzli.txt в следующем формате:

выделяем память под массив узлов(он и есть массив сетки)
считываем циклом значения узлов ниже от 1 до k+1
X1 X2 X3..Xk+1

-----
по идеее краевые условия можно хранить в этом же файле в следующем формате:
тип краевого условия на левой границе (1 - первые условия 2 - вторые)
тип краевого условия на правой границе (1 - первые условия 2 - вторые)

Краевых условий 3 типа у нас нет по условию!
/

int quantity_of_finite_elements = 0; // количество конечных элементов
int quantity_of_units = 0; // количество узлов

int type_of_left_boundary_con = 0; // тип левого краевого условия
int type_of_right_boundary_con = 0; // тип правого краевого условия
global_matrix_A = new double[quantity_of_units];

for (int i = 0; i < quantity_of_units; i++)
{
global_matrix_A = new double[quantity_of_units];
}

for (int i = 0; i < quantity_of_units; i++)
{
for (int j = 0; j < quantity_of_units; j++)
{
global_matrix_A[j] = 0;
}
}

vector_B = new double[quantity_of_units];
for (int i = 0; i < quantity_of_units; i++)
{
vector_B = 0;
}

grid = read_grid(grid, quantity_of_finite_elements, quantity_of_units);
read_boundary_cons(type_of_left_boundary_con, type_of_right_boundary_con);

global_compilator(quantity_of_units, quantity_of_finite_elements, grid);// генерация глобальной матрицы

global_matrix_A = chol_decomposition(quantity_of_units, global_matrix_A); // разложение Холецкого
switch (type_of_left_boundary_con) {
case 1:
set_left_boundary_con_1();
break;
case 2:
set_left_boundary_con_2();
break;
}
switch (type_of_right_boundary_con) {
case 1:
set_right_boundary_con_1(quantity_of_units);
break;
case 2:
set_right_boundary_con_2(quantity_of_units);
}
vector_B = ini_vect_y(quantity_of_units, vector_B, global_matrix_A); //прямой ход Гауса
vector_B = ini_vect_x(quantity_of_units, vector_B, global_matrix_A); // обратный ход Гауса

int s = 9;
}