#include "Solver.hpp"
void Solver::SolveImplicit(System& program, double tstop) const {
	std::ofstream output(filename_impl);
	//output << "t x y T" << std::endl;

	for (double t = 0.0; t < tstop; t += delta) {
		/* Обработка узлов по оси X */
		for (int i = 1; i < program.LineX().size() - 1; i++) {
			std::vector<Node*> temperature;
			Node* cur = program.LineX()[i];
			while (cur) {
				/* Проверка на существование узла */
				if (cur->r() && cur->r()->X() - cur->X() > program.step()) {
					temperature.push_back(cur);
					SolveLine(program, temperature);
					temperature.clear();
					cur = cur->r();
				}
				else {
					temperature.push_back(cur);
					cur = cur->r();
				}
			}
			SolveLine(program, temperature);
		}
		/* Обработка узлов по оси Y */
		for (int i = 1; i < program.LineY().size() - 1; i++) {
			std::vector<Node*> temperature;
			Node* cur = program.LineY()[i];
			while (cur) {
				/* Проверка на существование узла */
				if (cur->u() && cur->u()->Y() - cur->Y() > program.step()) {
					temperature.push_back(cur);
					SolveLine(program, temperature);
					temperature.clear();
					cur = cur->u();
				}
				else {
					temperature.push_back(cur);
					cur = cur->u();
				}
			}
			SolveLine(program, temperature);
		}
		for (auto line : program.Nodes()) {
			for (auto node : line)
				output << t+1 << ' ' << node->X() << ' ' << node->Y() << ' ' << node->T() << '\n';
		}
		output << "\n\n";
	}
}

void Solver::SolveExplicit(System& sys, double tstop) const {
	std::ofstream output(filename_expl);
	//output << "t x y T" << std::endl;

	for (double t = 0.; t < tstop; t += delta) {
		for (auto line : sys.Nodes())
			for (auto node : line) {
				/* Проверка на внутренний узел */
				if (!node->IsBound()) {
					/* Tx = T_right - 2T_current + T_left / delta_x ^ 2 */
					/* Ty = T_upper - 2T_current + T_down / delta_y ^ 2*/
					/* T_new = delta_t * a * (delta_x + delta_y) + T_current */

					double tx = (node->r()->T() - 2 * node->T() + node->l()->T()) / pow(sys.step(), 2);
					double ty = (node->u()->T() - 2 * node->T() + node->d()->T()) / pow(sys.step(), 2);
					double t = delta * (tx + ty) * sys.a() + node->T();
					node->SetT(t);
				}
			}
		for (auto line : sys.Nodes()) {
			for (auto node : line)
				output << t + 1  << ' ' << node->X() << ' ' << node->Y() << ' ' << node->T() << '\n';
		}
		output << "\n\n";
	}
}

void Solver::SolveLine(System& sys, std::vector<Node*>& n) const {
	int size = n.size() - 2;
	double mu1 = n.front()->Dist(n[1]) / sys.step();
	double mu2 = n.back()->Dist(n[n.size() - 2]) / sys.step();
	/* Защита от нуля */
	if (mu2 == 0.)
		mu2 = .1;
	double val2 = -(2 * sys.a()) / (pow(sys.step(), 2)) - 1 / delta;
	double val1 = sys.a() / (pow(sys.step(), 2));
	std::vector<std::vector<double>> _Temperature(size);
	std::vector<double> right(size);
	for (int i = 0; i < _Temperature.size(); i++)
		_Temperature[i].resize(3, 0.0);

	_Temperature[0][0] = -(2 * sys.a()) / (mu1 * pow(sys.step(), 2)) - 1 / delta; /* Первый узел по X */
	_Temperature[0][1] = (2 * sys.a()) / ((mu1 + 1) * pow(sys.step(), 2)); /* Первый узел по Y */
	_Temperature.back()[1] = (2 * sys.a()) / ((mu2 + 1) * pow(sys.step(), 2));
	_Temperature.back()[2] = -(2 * sys.a()) / (mu2 * pow(sys.step(), 2)) - 1 / delta;

	for (int i = 1; i < size - 1; i++) {
		_Temperature[i][0] = val1;
		_Temperature[i][1] = val2;
		_Temperature[i][2] = val1;
	}

	for (int i = 0; i < right.size(); i++)
		right[i] = -n[i + 1]->T() / delta;
	right.front() += -(2 * sys.a() * n.front()->T()) / (mu1 * (mu1 + 1) * pow(sys.step(), 2));
	right.back() += -(2 * sys.a() * n.back()->T()) / (mu2 * (mu2 + 1) * pow(sys.step(), 2));
	std::vector<double> tmp = ThomasMethod(_Temperature, right);
	for (int i = 0; i < tmp.size(); i++)
		n[i + 1]->SetT(tmp[i]);
}

/* Метод прогонки для численного решения СЛАУ */
std::vector<double> Solver::ThomasMethod(std::vector<std::vector<double>>& A, std::vector<double>& b) const {
	int row = b.size() - 1;

	std::vector<double>alpha(row);
	std::vector<double>beta(row);

	alpha[0] = -A[0][1] / A[0][0];
	beta[0] = b[0] / A[0][0];

	for (int i = 1; i < row; i++) {
		double a_i = A[i][0];
		double b_i = A[i][1];
		double c_i = A[i][2];
		double y_i = b_i + a_i * alpha[i - 1];
		
		alpha[i] = -c_i / y_i;
		beta[i] = (b[i] - a_i * beta[i - 1]) / y_i;
	}

	std::vector<double> result(b.size());

	result.back() = (b.back() - A.back()[1] * beta.back()) / (A.back()[2] + A.back()[1] * alpha.back());
	for (int i = row - 1; i > -1; i--)
		result[i] = alpha[i] * result[i + 1] + beta[i];
	return result;
}