solver.cpp

#include <iostream>
#include <iomanip>
#include <cmath>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <unordered_map>
#include "data.h"
#include "gauss.h"
#include "solver.h"

using namespace std;

typedef double number;

void test_main(GlobalData& globalData, vector<Node>& nodes, vector<Element>& elements, Grid& grid, vector<int>& borders, ElemUniv& elementUniv, int num_pc)
{
	initializeElementUniv(elementUniv, num_pc); //inicjalizacja elementu uniwersalnego

	//Obliczanie jakobianow oraz macierzy H
	cout << endl << "Obliczanie Jakobianow oraz macierzy H!" << endl;
	for (int i = 0; i < elements.size(); ++i)
	{
		//liczenie jakobianu dla kazdego elementu
		calculateJacobians(nodes, elements[i], elementUniv, num_pc);

		//liczenie macierzy H dla kazdego elementu
		calculate_matrix_H(elements[i], globalData, elementUniv, num_pc);
	}

	//Agregacja lokalnej macierzy H do macierzy globalnej H
	save_to_grid(grid, nodes, elements); //zapis do siatki
	cout << "Agregacja macierzy H do macierzy globalnej H -> ";
	aggregate_H(grid);
	cout << "H_Global(" << grid.H_Global.size() << "x" << grid.H_Global[0].size() << ")!" << endl;

	//Obliczanie macierzy Hbc i wektora P
	cout << endl << "Obliczanie macierzy Hbc i wektora P!" << endl;
	for (int i = 0; i < elements.size(); ++i)
	{
		//liczenie macierzy Hbc dla kazdego elementu
		calculate_matrix_Hbc(elements[i], globalData, nodes, elementUniv, num_pc);

		//zapisywanie Hbc do globanej H
		add_Hbc_to_H(elements[i]);

		//liczenie wektora P dla kazdego elementu
		calculate_P(elements[i], globalData, nodes, elementUniv, num_pc);
	}


	//Agregacja lokalnej macierzy H do macierzy globalnej H
	save_to_grid(grid, nodes, elements); //zapis do siatki!!!
	cout << "Agregacja wektora P do wektora globalnego P -> ";
	aggregate_P(grid);
	cout << "H_Global(" << grid.P_Global.size() << ")!" << endl;

	cout << fixed << setprecision(4);
	////SPRAWDZENIE WYNIKOW - JAKOBIAN
	// cout << endl << "Obliczanie Jakobianow" << endl;
	//cout << "Wyniki dla Elementu nr 1 (indeks 0):" << endl;
	//cout << "DetJ w punkcie calkowania 0: " << elements[0].jacobian[0].detJ << endl;
	//cout << "dN/dx w punkcie calkowania 0 (dla 4 wezlow): ";
	//for (int j = 0; j < 4; j++)
	//{
	//	cout << elements[0].jacobian[0].dN_dx[j] << " ";
	//}
	//cout << endl;
	//system("pause");

	////SPRAWDZANIE WYNIKOW - MACIERZ H
	//cout << endl << "Obliczanie macierzy H" << endl;
	//cout << "Wyniki dla Elementu nr 1 (indeks 0):" << endl;
	//for (int i = 0; i < npc; ++i)
	//{
	//	for (int j = 0; j < npc; ++j)
	//	{
	//		cout << setw(4) << elements[0].H[i][j] << " ";
	//	}
	//	cout << endl;
	//}
	//system("pause");

	////SPRAWDZENIE WYNIKOW - MACIERZ GLOBALNA H
	//cout << endl << "Macierz globalna H(" << grid.H_Global.size() << "x" << grid.H_Global[0].size() << ") : " << endl;
	//for (int i = 0; i < grid.H_Global.size(); ++i)
	//{
	//	for (int j = 0; j < grid.H_Global[i].size(); ++j)
	//	{
	//		cout << setw(8) << grid.H_Global[i][j] << " ";
	//	}
	//	cout << endl;
	//}
	//system("pause");

	////SPRAWDZANIE WYNIKOW - MACIERZ Hbc
	//cout << endl << "Obliczanie macierzy Hbc" << endl;
	//cout << "Wyniki dla Elementu nr 1 (indeks 0):" << endl;
	//for (int i = 0; i < npc; ++i)
	//{
	//	for (int j = 0; j < npc; ++j)
	//	{
	//		cout << setw(4) << elements[0].Hbc[i][j] << " ";
	//	}
	//	cout << endl;
	//}
	//system("pause");


	////SPRAWDZENIE WYNIKOW - WEKTOR GLOBALNY P
	//cout << endl << "Wektor globalny P(" << grid.P_Global.size() << ") : " << endl;
	//for (int i = 0; i < grid.P_Global.size(); ++i)
	//{
	//	cout << setw(8) << grid.P_Global[i] << " ";
	//}
	//cout << endl;
	//system("pause");

	return;
}

void Time_Calculation(GlobalData& globalData, vector<Node>& nodes, vector<Element>& elements, Grid& grid, vector<int>& borders, ElemUniv& elementUniv, int num_pc)
{

	initializeElementUniv(elementUniv, num_pc); //inicjalizacja elementu uniwersalnego
	cout << "Inicjalizacja elementu uniwersalnego zakonczona!" << endl << endl;

	//Obliczanie jakobianow oraz macierzy H
	cout << endl << "Obliczanie Jakobianow, macierzy H, macierzy Hbc, wektora P i macierzy C..." << endl << endl;
	for (int i = 0; i < elements.size(); ++i)
	{
		//liczenie jakobianu dla kazdego elementu
		calculateJacobians(nodes, elements[i], elementUniv, num_pc);
		//cout << "Obliczanie jakobianu dla elementu " << i + 1 << " zakonczone!" << endl;

		//liczenie macierzy H dla kazdego elementu
		calculate_matrix_H(elements[i], globalData, elementUniv, num_pc);
		//cout << "Obliczanie macierzy H dla elementu " << i + 1 << " zakonczone!" << endl;

		//liczenie macierzy Hbc dla kazdego elementu
		calculate_matrix_Hbc(elements[i], globalData, nodes, elementUniv, num_pc);
		//cout << "Obliczanie macierzy Hbc dla elementu " << i + 1 << " zakonczone!" << endl;

		//liczenie wektora P dla kazdego elementu
		calculate_P(elements[i], globalData, nodes, elementUniv, num_pc);
		//cout << "Obliczanie wektora P dla elementu " << i + 1 << " zakonczone!" << endl;

		//liczenie macierzy C dla kazdego elementu
		calculate_C(elements[i], globalData, elementUniv, num_pc);
		//cout << "Obliczanie macierzy C dla elementu " << i + 1 << " zakonczone!" << endl;

		//zapisywanie Hbc do globanej H
		add_Hbc_to_H(elements[i]);
	}
	cout << endl;

	save_to_grid(grid, nodes, elements); //zapis do siatki

	//Agregacja lokalnych macierzy (i wektora) do macierzy (i wektora) globalnych
	cout << "Agregacja lokalnych macierzy (i wektora) do macierzy (i wektora)..." << endl << endl;
	aggregate_H(grid);
	aggregate_P(grid);
	aggregate_C(grid);

	cout << "Obliczenia zakonczone!" << endl << endl << "Rozpoczynanie obliczania problemu niestacjonarnego..." << endl << endl;

	//Kroki czasowe
	vector<double> TempPrev(grid.nN.size(), globalData.InitialTemp); //wektor temperatur w nastepnym kroku czasowym
	vector<double> TempNext(grid.nN.size()); //wektor temperatur w nastepnym kroku czasowym

	double deltaT = globalData.TimeStep; //krok czasowy
	double totalTime = globalData.SimulationTime; //czas symulacji
	double currentTime = 0.0; //aktualny czas

	//Stan poczatkowy
	cout << "Time: " << currentTime << " s" << endl;
	cout << "Min: " << globalData.InitialTemp << " , Max: " << globalData.InitialTemp << endl << endl;

	while (currentTime < totalTime)
	{
		currentTime += deltaT;

		//A -> (H + C/delta_TimeStep)
		//B -> P + C/delta_TimeStep * Temp0
		//X -> Temperatury w nastepnym kroku czasowym (TempNext)

		vector<vector<double>> A(grid.nN.size(), vector<double>(grid.nN.size())); //macierz A
		vector<double> B(grid.nN.size()); //wektor B

		for (int i = 0; i < grid.nN.size(); ++i)
		{
			//Budowa wektora B
			//B[i] = P[i] + sum(C/deltaT * t0)
			double c_dt_t0 = 0.0; //co to?
			for (int j = 0; j < grid.nN.size(); ++j)
			{
				c_dt_t0 += (grid.C_Global[i][j] / deltaT) * TempPrev[j];
			}
			B[i] = grid.P_Global[i] + c_dt_t0;

			//Budowa macierzy A
			for (int j = 0; j < grid.nN.size(); ++j)
			{
				A[i][j] = grid.H_Global[i][j] + (grid.C_Global[i][j] / deltaT);
			}
		}

		if (currentTime == deltaT) // Wypisz tylko w 1. kroku
		{
			cout << "\n--- DIAGNOSTYKA PETLI (Czas: " << currentTime << ") ---" << endl;
			cout << "P_Global[0] wewnatrz petli: " << grid.P_Global[0] << endl;
			cout << "B[0] (Prawa Strona): " << B[0] << " (Oczekiwane: > 12000)" << endl;
			cout << "A[0][0] (Lewa Strona): " << A[0][0] << endl;
			cout << "TempPrev[0]: " << TempPrev[0] << endl;
			cout << "------------------------------------------\n";
		}

		solve_equation(A, B, TempNext, grid.nN.size()); //rozwiazanie ukladu rownan (plik solver.cpp)

		double minTemp = TempNext[0];
		double maxTemp = TempNext[0];
		for (double val : TempNext)
		{
			if (val < minTemp) minTemp = val;
			if (val > maxTemp) maxTemp = val;
		}

		cout << "Time: " << currentTime << " s" << endl;
		cout << "Min: " << minTemp << " , Max: " << maxTemp << endl << endl;

		TempPrev = TempNext; //aktualizacja temperatury poprzedniej
	}

	system("pause");

	return;
}