Добавьте файлы проекта.

src/Form.cpp

@@ -0,0 +1,43 @@
+#include "Form.hpp"
+
+size_t Form::counter_ = 0;
+
+double Form::Function(double, double) {
+    return 0;
+}
+
+std::pair<double, double> Form::Deriative(double, double) {
+    return { 0, 0 };
+}
+
+std::pair<double, double> Form::size() {
+    return { 0, 0 };
+}
+
+bool Form::Inhere(double, double) {
+    return false;
+}
+
+std::pair<double, double> Form::missX(double) {
+    return { 0, 0 };
+}
+
+std::pair<double, double> Form::missY(double) {
+    return { 0, 0 };
+}
+
+Form::Form() {
+    id_ = counter_++;
+    excluded_ = false;
+    _boundtype = -1;
+}
+
+bool Form::Excluded() const {
+    return excluded_;
+}
+
+int Form::GetB() { return _boundtype; }
+
+bool Form::operator==(size_t id) const {
+    return id_ == id;
+}

src/Mesh.cpp

@@ -0,0 +1,165 @@
+#include "Mesh.hpp"
+#include <iostream>
+#include <Object.hpp>
+
+Mesh::Mesh(Object& obj, double step) : _obj(obj), _step(step) {
+
+	for (double y = 0.0; y <= _obj.Height(); y += _step) {
+		_mesh.push_back(std::vector<Node*>());
+		for (double x = 0.0; x <= _obj.Width(); x += _step) {
+			_mesh.back().push_back(new Node(x, y));
+		}
+	}
+	LinkX();
+	LinkY();
+	Adapt();
+}
+
+void Mesh::LinkX() {
+
+	for (int i = 0; i < _mesh.size(); i++) {
+		_mesh[i][0]->LinkX(nullptr, _mesh[i][1]);
+		for (int j = 1; j < _mesh[i].size() - 1; j++)
+			_mesh[i][j]->LinkX(_mesh[i][j - 1], _mesh[i][j + 1]);
+		_mesh[i].back()->LinkX(_mesh[i][_mesh[i].size() - 2], nullptr);
+	}
+	for (int i = 0; i < _mesh.size(); i++)
+		_hlines.push_back(_mesh[i][0]);
+}
+
+
+void Mesh::LinkY() {
+	for (int j = 0; j < _mesh[0].size(); j++) {
+		_mesh[0][j]->LinkY(nullptr, _mesh[1][j]);
+		for (int i = 1; i < _mesh.size() - 1; i++)
+			_mesh[i][j]->LinkY(_mesh[i - 1][j], _mesh[i + 1][j]);
+		_mesh[_mesh.size() - 1][j]->LinkY(_mesh[_mesh.size() - 2][j], nullptr);
+	}
+	for (int i = 0; i < _mesh[0].size(); i++)
+		_vlines.push_back(_mesh[0][i]);
+}
+
+
+void Mesh::Adapt() {
+	for (int i = 0; i < _mesh.size(); i++) {
+		int s = _mesh[i].size();
+		for (int j = 0; j < s; j++) {
+			if (!_obj.Inhere(_mesh[i][j]->X(), _mesh[i][j]->Y())) {
+				Delnode(i, j);
+				j--;
+				s--;
+			}
+		}
+	}
+}
+
+void Mesh::ShowLinks() {
+	for (auto line : _mesh) {
+		for (auto node : line) {
+			if (node->d())
+				std::cout << "|  ";
+		}
+		std::cout << '\n';
+		for (auto node : line) {
+			if (node->l()) {
+				std::cout << '-';
+			}
+			std::cout << 'N';
+			if (node->r()) {
+				std::cout << '-';
+			}
+			else {
+				std::cout << '\n';
+			}
+		}
+		for (auto node : line) {
+			if (node->u())
+				std::cout << "|";
+			std::cout << "  ";
+
+		}
+		std::cout << '\n';
+	}
+}
+
+void Mesh::Delnode(int i, int j) {
+	Node* node = _mesh[i][j];
+	double bndX1 = _obj.Fillx(node->X(), node->Y()).first;
+	double bndX2 = _obj.Fillx(node->X(), node->Y()).second;
+	double bndY1 = _obj.Filly(node->X(), node->Y()).first;
+	double bndY2 = _obj.Filly(node->X(), node->Y()).second;
+	int btype = _obj.Who(node->X(), node->Y())->GetB();
+	if (node->l()) {
+		if (node->l()->X() != bndX2 && node->l()->X() != bndX1) {
+			if (bndX1 != bndX2) {
+				Node* left = new Node(bndX2, node->Y(), btype);
+				Node* right = new Node(bndX1, node->Y(), btype);
+				node->l()->r() = left;
+				if (node->r())
+					node->r()->l() = right;
+				left->LinkX(node->l(), right);
+				right->LinkX(left, node->r());
+				node->l() = right;
+				_mesh[i].push_back(left);
+				_mesh[i].push_back(right);
+			}
+			else {
+				Node* left = new Node(bndX2, node->Y(), btype);
+				node->l()->r() = left;
+				if (node->r())
+					node->r()->l() = left;
+				left->LinkX(node->l(), node->r());
+				node->l() = left;
+				_mesh[i].push_back(left);
+			}
+		}
+		else
+			node->l()->r() = node->r();
+	}
+	if (node->r()) {
+		node->r()->l() = node->l();
+	}
+	if (node->d()) {
+		if (node->d()->Y() != bndY2 && node->d()->Y() != bndY1) {
+			if (bndY2 != bndY1) {
+				Node* down = new Node(node->X(), bndY2, btype);
+				Node* up = new Node(node->X(), bndY1, btype);
+				node->d()->u() = down;
+				if (node->u())
+					node->u()->d() = up;
+				down->LinkY(node->d(), up);
+				up->LinkY(down, node->u());
+				node->d() = up;
+				_mesh[i].push_back(down);
+				_mesh[i].push_back(up);
+			}
+			else {
+				Node* down = new Node(node->X(), bndY2, btype);
+				node->d()->u() = down;
+				if (node->u())
+					node->u()->d() = down;
+				down->LinkY(node->d(), node->u());
+				node->d() = down;
+				_mesh[i].push_back(down);
+			}
+		}
+		else
+			node->d()->u() = node->u();
+	}
+	if (node->u()) {
+		node->u()->d() = node->d();
+	}
+	_mesh[i].erase(_mesh[i].begin() + j);
+	delete node;
+}
+std::vector<std::vector<Node*>>& Mesh::Nodes() { return _mesh; }
+std::vector<Node*>& Mesh::LineX() { return _hlines; }
+std::vector<Node*>& Mesh::LineY() { return _vlines; }
+Mesh::~Mesh() {
+	for (auto line : _mesh)
+		for (auto node : line)
+			delete node;
+}
+
+
+

src/Node.cpp

@@ -0,0 +1,87 @@
+#include "Node.hpp"
+#include <iostream>
+double Node::X() const { return _x; }
+double Node::Y() const { return _y; }
+double Node::T() const {
+	if (_btype == 1)
+		return 100.;
+
+	if (_btype == 2) {
+		if (!_left)
+			if (_right)
+				return _right->T();
+		if (!_right)
+			if (_left)
+				return _left->T();
+		if (!_above)
+			if (_below)
+				return _below->T();
+		if (!_below)
+			if (_above)
+				return _above->T();
+		if (_right && _left) {
+			if (_right->IsBound())
+				return _left->T();
+			return _right->T();
+		}
+		if (_above && _below) {
+			if (_above->IsBound())
+				return _below->T();
+			return _above->T();
+		}
+	}
+
+	if (_btype == 3) {
+		if (!_left)
+			if (_right)
+				return _right->T() / (1 + Dist(_right));
+		if (!_right)
+			if (_left)
+				return _left->T() / (1 + Dist(_left));
+		if (!_above)
+			if (_below)
+				return _below->T() / (1 + Dist(_below));
+		if (!_below)
+			if (_above)
+				return _above->T() / (1 + Dist(_above));
+		if (_right && _left) {
+			if (_right->IsBound())
+				return _left->T() / (1 + Dist(_left));
+			return _right->T() / (1 + Dist(_right));
+		}
+		if (_above && _below) {
+			if (_above->IsBound())
+				return _below->T() / (1 + Dist(_below));
+			return _above->T() / (1 + Dist(_above));
+		}
+	}
+
+	return _t;
+}
+
+double Node::Dist(const Node* to) const {
+	return std::sqrt(pow(X() - to->X(), 2) + pow(Y() - to->Y(), 2));
+}
+
+Node*& Node::l() { return _left; }
+Node*& Node::r() { return _right; }
+Node*& Node::u() { return _above; }
+Node*& Node::d() { return _below; }
+
+void Node::LinkX(Node* l, Node* r) {
+	_left = l;
+	_right = r;
+}
+
+void Node::LinkY(Node* d, Node* u) {
+	_below = d;
+	_above = u;
+}
+
+void Node::SetT(double t) {
+	_t = t;
+}
+
+bool Node::IsBound() { return _btype; }
+void Node::SetB(int type) { _btype = type; }
+

src/Object.cpp

@@ -0,0 +1,97 @@
+#include "Object.hpp"
+
+Object::Object() : _w(0), _h(0) {}
+
+double Object::Inhere(double x, double y) {
+    for (auto form : forms_) {
+        if (form->Excluded()) {
+            if (form->Inhere(x, y)) {
+                return false;
+            }
+        }
+        else {
+            if (form->Inhere(x, y)) {
+                return true;
+            }
+        }
+    }
+    return false;
+}
+
+void Object::Updsize() {
+    for (auto form : forms_) {
+        _w = std::max(_w, form->size().first);
+        _h = std::max(_h, form->size().second);
+    }
+}
+
+bool Object::Add_Form(const std::string& name, std::map<std::string, double>& args, bool excluded, int btype) {
+    if (name == "Rectangle") {
+        forms_.push_back(new Rectangle(args["a"], args["b"], args["h_x"], args["h_y"], excluded, btype));
+        Updsize();
+        return true;
+    }
+    else if (name == "Circle") {
+        forms_.push_back(new Circle(args["a"], args["b"], args["h_x"], args["h_y"], excluded, btype));
+        Updsize();
+        return true;
+    }
+    else if (name == "Arc") {
+        forms_.push_back(new Arc(args["a"], args["b"], args["h_x"], args["h_y"], excluded, btype));
+        Updsize();
+        return true;
+    }
+    return false;
+}
+
+
+
+std::pair<double, double> Object::Fillx(double x, double y) {
+    for (auto form : forms_) {
+        if (form->Inhere(x, y)) {
+            return form->missX(y);
+        }
+    }
+    return { 0, 0 };
+}
+
+std::pair<double, double> Object::Filly(double x, double y) {
+    for (auto form : forms_) {
+        if (form->Inhere(x, y)) {
+            return form->missY(x);
+        }
+    }
+    return { 0, 0 };
+}
+
+double Object::Width() const {
+    return _w;
+}
+
+double Object::Height() const {
+    return _h;
+}
+
+//std::vector<size_t> Object::Get_IDs() {
+//    std::vector<size_t> ids;
+//    ids.reserve(forms_.size());
+//    for (auto form : forms_) {
+//        ids.push_back(form->Get_ID());
+//    }
+//    return ids;
+//}
+
+Form* Object::Who(double x, double y) {
+    for (auto form : forms_) {
+        if (form->Inhere(x, y)) {
+            return form;
+        }
+
+    }
+    return forms_.back();
+}
+Object::~Object() {
+    for (auto form : forms_)
+        delete form;
+}
+

src/Primitives.cpp

@@ -0,0 +1,98 @@
+#include "Primitives.hpp"
+
+Rectangle::Rectangle(double a, double b, double h_x, double h_y, bool excluded, int btype) : a_(a), b_(b), h_x_(h_x), h_y_(h_y) {
+    excluded_ = excluded;
+    _boundtype = btype;
+}
+
+std::pair<double, double> Rectangle::missX(double y) {
+    return { 0.5 / h_x_ + a_, -0.5 / h_x_ + a_ };
+}
+std::pair<double, double> Rectangle::missY(double x) {
+    return { 0.5 / h_y_ + b_, -0.5 / h_y_ + b_ };
+}
+
+std::pair<double, double> Rectangle::size() {
+    return { 1 / h_x_, 1 / h_y_ };
+}
+
+double Rectangle::Function(double x, double y) {
+    return std::max(h_x_ * std::abs(x - a_), h_y_ * std::abs(y - b_));
+}
+
+std::pair<double, double> Rectangle::Deriative(double x, double y) {
+    return { (h_x_ / 2) * ((x - a_) / std::abs(x - a_)), (h_y_ / 2) * ((y - b_) / std::abs(y - b_)) };
+}
+
+bool Rectangle::Inhere(double x, double y) {
+    return Function(x, y) <= 0.5;
+}
+
+Circle::Circle(double a, double b, double h_x, double h_y, bool excluded, int btype) : a_(a), b_(b), h_x_(h_x), h_y_(h_y) {
+    excluded_ = excluded;
+    _boundtype = btype;
+}
+
+std::pair<double, double> Circle::missY(double x) {
+    return { std::sqrt(1 - pow((h_x_ * (x - a_)), 2)) / h_y_ + b_, -std::sqrt(1 - pow((h_x_ * (x - a_)), 2)) / h_y_ + b_ };
+}
+std::pair<double, double> Circle::missX(double y) {
+    return { std::sqrt(1 - pow((h_y_ * (y - b_)), 2)) / h_x_ + a_, -std::sqrt(1 - pow((h_y_ * (y - b_)), 2)) / h_x_ + a_ };
+}
+
+double Circle::Function(double x, double y) {
+    return pow(h_x_ * (x - a_), 2) + pow(h_y_ * (y - b_), 2);
+}
+
+std::pair<double, double> Circle::Deriative(double x, double y) {
+    return { 2 * h_x_ * (x - a_), 2 * h_y_ * (y - b_) };
+}
+
+std::pair<double, double> Circle::size() {
+    return { 1 / h_x_, 1 / h_y_ };
+}
+
+bool Circle::Inhere(double x, double y) {
+    return Function(x, y) <= 1;
+}
+
+Arc::Arc(double a, double b, double h_x, double h_y, bool excluded, int btype) : a_(a), b_(b), h_x_(h_x), h_y_(h_y) {
+    excluded_ = excluded;
+    _boundtype = btype;
+}
+
+std::pair<double, double> Arc::missY(double x) {
+    return { std::sqrt(1 - pow((h_x_ * (x - a_)), 2)) / h_y_ + b_, std::sqrt(1 - pow((h_x_ * (x - a_)), 2)) / h_y_ + b_ };
+}
+std::pair<double, double> Arc::missX(double y) {
+    return { std::sqrt(1 - pow((h_y_ * (y - b_)), 2)) / h_x_ + a_, std::sqrt(1 - pow((h_y_ * (y - b_)), 2)) / h_x_ + a_ };
+}
+
+double Arc::Function(double x, double y) {
+    if (x >= a_ && y >= b_) {
+        return pow(h_x_ * (x - a_), 2) + pow(h_y_ * (y - b_), 2);
+    }
+    return -1.0;
+}
+
+std::pair<double, double> Arc::Deriative(double x, double y) {
+    if (x >= a_ && y >= b_) {
+        return { 2 * h_x_ * (x - a_), 2 * h_y_ * (y - b_) };
+    }
+    if (x < a_) {
+        std::cout << "x < a\n";
+    }
+    if (y < b_) {
+        std::cout << "y < b\n";
+    }
+    return { -1.0, -1.0 };
+}
+
+std::pair<double, double> Arc::size() {
+    return { 1 / h_x_, 1 / h_y_ };
+}
+
+bool Arc::Inhere(double x, double y) {
+    return Function(x, y) >= 1;
+}
+

src/Solver.cpp

@@ -0,0 +1,133 @@
+#include "Solver.hpp"
+void Solver::SolveExplicit(System& program, double tstop) const {
+	std::ofstream ExplicitOut(_name_1);
+	for (double t = 0.0; t < tstop; t += delta) {
+		/* <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD> <20><> <20><><EFBFBD> X */
+		for (int i = 1; i < program.LineX().size() - 1; i++) {
+			std::vector<Node*> temperature;
+			Node* cur = program.LineX()[i];
+			while (cur) {
+				/* <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD> */
+				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);
+		}
+		/* <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD> <20><> <20><><EFBFBD> Y */
+		for (int i = 1; i < program.LineY().size() - 1; i++) {
+			std::vector<Node*> temperature;
+			Node* cur = program.LineY()[i];
+			while (cur) {
+				/* <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD> */
+				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)
+				ExplicitOut << node->X() << ' ' << node->Y() << ' ' << node->T() << '\n';
+		}
+		ExplicitOut << "\n\n";
+	}
+}
+
+void Solver::SolveImplicit(System& sys, double tstop) const {
+	std::ofstream EmplicitOut(_name_2);
+	for (double t = 0.; t < tstop; t += delta) {
+		for (auto line : sys.Nodes())
+			for (auto node : line) {
+				/* <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD> */
+				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 
+					(<28><><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> a = 1) */
+
+					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) + node->T();
+					node->SetT(t);
+				}
+			}
+		for (auto line : sys.Nodes()) {
+			for (auto node : line)
+				EmplicitOut << node->X() << ' ' << node->Y() << ' ' << node->T() << '\n';
+		}
+		EmplicitOut << "\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();
+	/* <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><> <20><><EFBFBD><EFBFBD> */
+	if (mu2 == 0.)
+		mu2 = .1;
+	double val2 = -(2 * sys.a1()) / (pow(sys.step(), 2)) - 1 / delta;
+	double val1 = sys.a1() / (pow(sys.step(), 2));
+	std::vector<std::vector<double>> next(size);
+	std::vector<double> right(size);
+	for (int i = 0; i < next.size(); i++)
+		next[i].resize(3, 0.0);
+	next[0][0] = -(2 * sys.a1()) / (mu1 * pow(sys.step(), 2)) - 1 / delta;//val2;
+	next[0][1] = (2 * sys.a1()) / ((mu1 + 1) * pow(sys.step(), 2));// val1;
+	next.back()[1] = (2 * sys.a1()) / ((mu2 + 1) * pow(sys.step(), 2));// val1;
+	next.back()[2] = -(2 * sys.a1()) / (mu2 * pow(sys.step(), 2)) - 1 / delta;//val2;
+	for (int i = 1; i < size - 1; i++) {
+		next[i][0] = val1;
+		next[i][1] = val2;
+		next[i][2] = val1;
+	}
+
+	for (int i = 0; i < right.size(); i++)
+		right[i] = -n[i + 1]->T() / delta;
+	right.front() += -(2 * sys.a1() * n.front()->T()) / (mu1 * (mu1 + 1) * pow(sys.step(), 2));
+	right.back() += -(2 * sys.a1() * n.back()->T()) / (mu2 * (mu2 + 1) * pow(sys.step(), 2));
+	std::vector<double> tmps = ThomasMethod(next, right);
+	for (int i = 0; i < tmps.size(); i++)
+		n[i + 1]->SetT(tmps[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 = A[i][0];
+		double b1 = A[i][1];
+		double c = A[i][2];
+		alpha[i] = -c / (a * alpha[i - 1] + b1);
+		beta[i] = (b[i] - a * beta[i - 1]) / (a * alpha[i - 1] + b1);
+	}
+
+	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;
+}

src/System.cpp

@@ -0,0 +1,30 @@
+#include "System.hpp"
+
+void System::DefineBounds(int l, int t, int r, int b) {
+	Node* cur = _mesh.LineX().front();
+	while (cur) {
+		cur->SetB(b);
+		cur = cur->r();
+	}
+	cur = _mesh.LineX().back();
+	while (cur) {
+		cur->SetB(t);
+		cur = cur->r();
+	}
+	cur = _mesh.LineY().front();
+	while (cur) {
+		cur->SetB(l);
+		cur = cur->u();
+	}
+	cur = _mesh.LineY().back();
+	while (cur->u()) {
+		cur->SetB(r);
+		cur = cur->u();
+	}
+}
+
+std::vector<std::vector<Node*>>& System::Nodes() { return _mesh.Nodes(); }
+std::vector<Node*>& System::LineX() { return _mesh.LineX(); }
+std::vector<Node*>& System::LineY() { return _mesh.LineY(); }
+
+

src/main.cpp

@@ -0,0 +1,99 @@
+#include <iostream>
+#include <fstream>
+#include <Solver.hpp>
+
+#define WIDTH 500.
+#define HEIGHT 400.
+
+#define ARC_RADIUS 150.
+
+//#define SQUARE_X 355.
+//#define SQUARE_Y 155.
+//#define SQUARE_SIDE 100.
+
+#define HOLE_X 355.
+#define HOLE_Y 155.
+#define HOLE_RADIUS 50.
+
+void visualize(std::ofstream& file, std::string filename, int time_end) {
+	file << "set cbrange [" << 0 << ":" << 100 << "]" << std::endl;
+	file << "set size ratio " << float(400) / 500 << "\nunset key\n" << "\nset palette defined (0 0 0 1, 0.25 0 1 1, 0.5 0 1 0, 0.75 1 1 0, 1 1 0 0)\n" << std::endl;
+	file << "do for [i=0:" << time_end - 1 << "]{" << std::endl;
+	file << "plot '" << filename << "' u 1:2:3 index i w points pt 5 palette" << std::endl;
+	file << "pause " << 0.000000001 << "}" << std::endl;
+	file << "pause mouse";
+}
+
+int main()
+{
+	int l = 3; //Bound conditions
+	int t = 3;
+	int r = 1;
+	int b = 3;
+	int r2 = 4;
+	int s = 3;
+
+	double step = 10; // Mesh step
+	double step2 = 5; // Mesh step
+
+	double time_step = 1;
+	double time_end = 100;
+
+	double C = 50.; // Material props
+
+	std::map<std::string, double> plate{ 
+		{"a", WIDTH / 2}, {"b", HEIGHT / 2}, {"h_x", 1 / WIDTH}, {"h_y", 1 / HEIGHT} 
+	};
+	std::map<std::string, double> arc{
+		{"a", WIDTH - ARC_RADIUS}, {"b", HEIGHT - ARC_RADIUS}, {"h_x", 1 / ARC_RADIUS}, {"h_y", 1 / ARC_RADIUS} 
+	};
+	std::map<std::string, double> circle{
+		{"a", HOLE_X}, {"b", HOLE_Y}, {"h_x", 1 / HOLE_RADIUS}, {"h_y", 1 / HOLE_RADIUS}
+	};
+
+	Object obj;
+	obj.Add_Form("Circle", circle, true, s);
+	obj.Add_Form("Arc", arc, true, r2);
+	obj.Add_Form("Rectangle", plate, false, 1);
+
+	System explicit10(obj, step, C);
+	explicit10.DefineBounds(l, t, r, b);
+
+	System explicit5(obj, step2, C);
+	explicit5.DefineBounds(l, t, r, b);
+
+	System implicit10(obj, step, C);
+	implicit10.DefineBounds(l, t, r, b);
+
+	System implicit5(obj, step2, C);
+	implicit5.DefineBounds(l, t, r, b);
+
+
+	Solver slv10("explicit10.txt", "implicit10.txt", time_step);
+	slv10.SolveExplicit(explicit10, time_end);
+	slv10.SolveImplicit(implicit10, time_end);
+
+	Solver slv5("explicit5.txt", "implicit5.txt", time_step);
+	slv5.SolveExplicit(explicit5, time_end);
+	slv5.SolveImplicit(implicit5, time_end);
+	std::ofstream script("es10.plt");
+	visualize(script, "explicit10.txt", time_end);
+	script.close();
+
+	script.open("is10.plt");
+	visualize(script, "implicit10.txt", time_end);
+	script.close();
+
+	script.open("es5.plt");
+	visualize(script, "explicit5.txt", time_end);
+	script.close();
+
+	script.open("is5.plt");
+	visualize(script, "implicit5.txt", time_end);
+	script.close();
+	return 0;
+}
+
+//std::map<std::string, double> square{ 
+	//	{"a", SQUARE_X}, {"b", SQUARE_Y}, {"h_x", 1 / SQUARE_SIDE}, {"h_y", 1 / SQUARE_SIDE} 
+	//};