gctl/lib/io/triangle_io.h

/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#ifndef _GCTL_TRIANGLE_IO_H
#define _GCTL_TRIANGLE_IO_H

// library's head file
#include "../core.h"
#include "../utility.h"
#include "../geometry.h"

namespace gctl
{
	/**
	 * @brief      Read Triangle's .node file into a 2D vertex array.
	 * （从 Triangle 程序的 .node 文件读入二维顶点数组）
	 * 
	 * First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)>
	 * Remaining lines: <vertex #> <x> <y> [attributes] [boundary marker]
	 * Blank lines and comments prefixed by `#' may be placed anywhere. Vertices must be numbered consecutively, starting from one or zero.
	 * The attributes, which are typically floating-point values of physical quantities 
	 * (such as mass or conductivity) associated with the nodes of a finite element mesh, 
	 * are copied unchanged to the output mesh. If -q, -a, -u, or -s is selected, 
	 * each new Steiner point added to the mesh will have quantities assigned to it by linear interpolation.
	 * If the fourth entry of the first line is `1', the last column of the remainder of the file is 
	 * assumed to contain boundary markers. Boundary markers are used to identify boundary vertices 
	 * and vertices resting on PSLG segments. The .node files produced by Triangle contain boundary 
	 * markers in the last column unless they are suppressed by the -B switch.
	 * 
	 * @param[in]  filename     Input filename without any extensions.（文件名，无后缀）
	 * @param      out_node     Output 2D vertex array. Returned by quote.（返回的二维顶点数组）
	 * @param[in]  packed       Indicates whether the starting index of edges is zero or not. The 
	 * edge's ordering is assumed to be start with one if this option is set to NotPacked. 
	 * The default value is Packed.（.node 文件中顶点序号从 0 开始则为真，为 1 则为假）
	 * @param      marker       If the pointer is not null. Return an int array of vertexs' marker.（指针不为空时返回一个包含顶点标记的整形数组）
	 * @param      attri        If the pointer is not null. Return a 2D integer array that contains
	 * attributes of the vertice.（指针不为空时返回一个二维数组，表示对应点的属性值）
	 */
	template <typename A>
	void read_Triangle_node(std::string filename, array<vertex<point2dc, A>> &out_node, index_packed_e packed = Packed, 
		array<int> *marker = nullptr, _2d_matrix *attri = nullptr)
	{
		std::ifstream tetin;
		open_infile(tetin, filename, ".node");

		std::string err_str, tmp_str;
		std::stringstream tmp_ss;
		int node_size = 0, node_dimen = 2, attri_num = 0, boundary_mark = 0;

		// read head info
		do
		{
			std::getline(tetin, tmp_str);
			if (tmp_str[0] != '#')
			{
				str2ss(tmp_str, tmp_ss);
				tmp_ss >> node_size >> node_dimen >> attri_num >> boundary_mark;
				if (tmp_ss.fail() || node_dimen != 2 || 
					(boundary_mark != 0 && boundary_mark != 1) || 
					attri_num < 0)
				{
					err_str = "Wrong head-info found in " + filename + 
						".node. From void gctl::read_Triangle_node(...)";
					throw runtime_error(err_str);
				}
			}
		}
		while (tmp_str[0] == '#');

		if (attri_num != 0 && attri != nullptr)
		{
			attri->resize(node_size, attri_num);
		}
		//else if (attri_num == 0 && attri != nullptr)
		//{
		//	err_str = "No attributes found in " + filename + 
		//		".node. From void gctl::read_Triangle_node(...)";
		//	throw runtime_error(err_str);
		//}

		if (boundary_mark != 0 && marker != nullptr)
		{
			marker->resize(node_size);
		}
		//else if (boundary_mark == 0 && marker != nullptr)
		//{
		//	err_str = "No markers found in " + filename + 
		//		".node. From void gctl::read_Triangle_node(...)";
		//	throw runtime_error(err_str);
		//}

		vertex<point2dc, A> tmp_vert;
		out_node.resize(node_size);
		if (boundary_mark+attri_num)
		{
			array<double> attri_line(boundary_mark+attri_num);
			_2d_matrix tmp_attri(node_size, boundary_mark+attri_num);
			while (std::getline(tetin, tmp_str))
			{
				if (tmp_str[0] != '#')
				{
					str2ss(tmp_str, tmp_ss);
					tmp_ss >> tmp_vert.id >> tmp_vert.x >> tmp_vert.y;
					for (int a = 0; a < boundary_mark+attri_num; a++)
						tmp_ss >> attri_line[a];

					if (tmp_ss.fail())
					{
						err_str = "Wrong file format found in " + filename + 
							".node. From void gctl::read_Triangle_node(...)";
						throw runtime_error(err_str);
					}

					if (packed == NotPacked) tmp_vert.id -= 1;
					out_node[tmp_vert.id] = tmp_vert;
					for (int a = 0; a < boundary_mark+attri_num; a++)
						tmp_attri[tmp_vert.id][a] = attri_line[a];
				}
			}

			if (attri != nullptr && attri_num)
			{
				for (int i = 0; i < node_size; i++)
				{
					for (int a = 0; a < attri_num; a++)
					{
						attri->at(i, a) = tmp_attri[i][a];
					}
				}
			}

			if (marker != nullptr && boundary_mark)
			{
				for (int i = 0; i < node_size; i++)
				{
					marker->at(i) = (int) tmp_attri[i][attri_num];
				}
			}
		}
		else
		{
			while (std::getline(tetin, tmp_str))
			{
				if (tmp_str[0] != '#')
				{
					str2ss(tmp_str, tmp_ss);
					tmp_ss >> tmp_vert.id >> tmp_vert.x >> tmp_vert.y;

					if (tmp_ss.fail())
					{
						err_str = "Wrong file format found in " + filename + 
							".node. From void gctl::read_Triangle_node(...)";
						throw runtime_error(err_str);
					}

					if (packed == NotPacked) tmp_vert.id -= 1;
					out_node[tmp_vert.id] = tmp_vert;
				}
			}
		}

		tetin.close();
		return;
	}

	/**
	 * @brief      Read Triangle's .ele file into a 2D triangle array.
	 * 
	 * First line: <# of triangles> <nodes per triangle> <# of attributes>
	 * Remaining lines: <triangle #> <node> <node> <node> ... [attributes]
	 * Blank lines and comments prefixed by `#' may be placed anywhere. Triangles must be numbered 
	 * consecutively, starting from one or zero. Nodes are indices into the corresponding .node file. 
	 * The first three nodes are the corner vertices, and are listed in counterclockwise order 
	 * around each triangle. (The remaining nodes, if any, depend on the type of finite element used.)
	 * As in .node files, the attributes are typically floating-point values of physical quantities 
	 * (such as mass or conductivity) associated with the elements (triangles) of a finite element mesh. 
	 * Because there is no simple mapping from input to output triangles, an attempt is made to 
	 * interpolate attributes, which may result in a good deal of diffusion of attributes among nearby 
	 * triangles as the triangulation is refined. Attributes do not diffuse across segments, 
	 * so attributes used to identify segment-bounded regions remain intact.
	 * In output .ele files, all triangles have three nodes each unless the -o2 switch is used, in 
	 * which case subparametric quadratic elements with six nodes are generated. The fourth, fifth, 
	 * and sixth nodes lie on the midpoints of the edges opposite the first, second, and third vertices.
	 *
	 * @param      out_tri   Output 2D triangle array. Returned by quote.
	 * @param[in]  filename  Input filename without any extensions.
	 * @param[in]  in_node   2D vertex array which contains vertex's information of the triangles.
	 * @param[in]  packed    Indicates whether the starting index of edges is zero or not. The 
	 * edge's ordering is assumed to be start with one if this option is set to NotPacked. 
	 * The default value is Packed.
	 * @param      attri     If the pointer is not null. Return a 2D integer array that contains
	 * attributes of the vertice.
	 */
	template <typename E, typename A>
	void read_Triangle_element(std::string filename, array<type_triangle2d<E>> &out_tri, 
		const array<vertex<point2dc, A>> &in_node, index_packed_e packed = Packed, 
		_2d_matrix *attri = nullptr)
	{
		std::ifstream tetin;
		open_infile(tetin, filename, ".ele");

		std::string err_str, tmp_str;
		std::stringstream tmp_ss;
		int ele_size = 0, node_per_ele = 0, attri_num = 0;

		// read head info
		do
		{
			std::getline(tetin, tmp_str);
			if (tmp_str[0] != '#')
			{
				str2ss(tmp_str, tmp_ss);
				tmp_ss >> ele_size >> node_per_ele >> attri_num;
				if (tmp_ss.fail() || node_per_ele != 3 || attri_num < 0)
				{
					err_str = "Wrong head-info found in " + filename + 
						".ele. From void gctl::read_Triangle_element(...)";
					throw runtime_error(err_str);
				}
			}
		}
		while (tmp_str[0] == '#');

		if (attri_num != 0 && attri != nullptr)
		{
			attri->resize(ele_size, attri_num);
		}
		//else if (attri_num == 0 && attri != nullptr)
		//{
		//	err_str = "No attributes found in " + filename + 
		//		".ele. From void gctl::read_Triangle_element(...)";
		//	throw runtime_error(err_str);
		//}

		int tmp_int, tmp_order[3];
		out_tri.resize(ele_size);
		if (attri_num)
		{
			array<double> attri_line(attri_num);
			_2d_matrix tmp_attri(ele_size, attri_num);
			while (std::getline(tetin, tmp_str))
			{
				if (tmp_str[0] != '#')
				{
					str2ss(tmp_str, tmp_ss);
					tmp_ss >> tmp_int >> tmp_order[0] >> tmp_order[1] >> tmp_order[2];
					for (int a = 0; a < attri_num; a++)
						tmp_ss >> attri_line[a];

					if (tmp_ss.fail())
					{
						err_str = "Wrong file format found in " + filename + 
							".ele. From void gctl::read_Triangle_element(...)";
						throw runtime_error(err_str);
					}

					if (packed == NotPacked)
					{
						tmp_int -= 1;
						out_tri[tmp_int].id = tmp_int;
						for (int j = 0; j < 3; j++)
						{
							out_tri[tmp_int].vert[j] = in_node.get(tmp_order[j]-1);
						}
					}
					else
					{
						out_tri[tmp_int].id = tmp_int;
						for (int j = 0; j < 3; j++)
						{
							out_tri[tmp_int].vert[j] = in_node.get(tmp_order[j]);
						}
					}

					for (int a = 0; a < attri_num; a++)
						tmp_attri[tmp_int][a] = attri_line[a];
				}
			}

			if (attri != nullptr)
			{
				for (int i = 0; i < ele_size; i++)
				{
					for (int a = 0; a < attri_num; a++)
					{
						attri->at(i, a) = tmp_attri[i][a];
					}
				}
			}
		}
		else
		{
			while (std::getline(tetin, tmp_str))
			{
				if (tmp_str[0] != '#')
				{
					str2ss(tmp_str, tmp_ss);
					tmp_ss >> tmp_int >> tmp_order[0] >> tmp_order[1] >> tmp_order[2];

					if (tmp_ss.fail())
					{
						err_str = "Wrong file format found in " + filename + 
							".ele. From void gctl::read_Triangle_element(...)";
						throw runtime_error(err_str);
					}

					if (packed == NotPacked)
					{
						tmp_int -= 1;
						out_tri[tmp_int].id = tmp_int;
						for (int j = 0; j < 3; j++)
						{
							out_tri[tmp_int].vert[j] = in_node.get(tmp_order[j]-1);
						}
					}
					else
					{
						out_tri[tmp_int].id = tmp_int;
						for (int j = 0; j < 3; j++)
						{
							out_tri[tmp_int].vert[j] = in_node.get(tmp_order[j]);
						}
					}
				}
			}
		}

		tetin.close();
		return;
	}

	/**
	 * @brief      Read Triangle's .ele file into a second order 2D triangle array.
	 * 
	 * First line: <# of triangles> <nodes per triangle> <# of attributes>
	 * Remaining lines: <triangle #> <node> <node> <node> ... [attributes]
	 * Blank lines and comments prefixed by `#' may be placed anywhere. Triangles must be numbered 
	 * consecutively, starting from one or zero. Nodes are indices into the corresponding .node file. 
	 * The first three nodes are the corner vertices, and are listed in counterclockwise order 
	 * around each triangle. (The remaining nodes, if any, depend on the type of finite element used.)
	 * As in .node files, the attributes are typically floating-point values of physical quantities 
	 * (such as mass or conductivity) associated with the elements (triangles) of a finite element mesh. 
	 * Because there is no simple mapping from input to output triangles, an attempt is made to 
	 * interpolate attributes, which may result in a good deal of diffusion of attributes among nearby 
	 * triangles as the triangulation is refined. Attributes do not diffuse across segments, 
	 * so attributes used to identify segment-bounded regions remain intact.
	 * In output .ele files, all triangles have three nodes each unless the -o2 switch is used, in 
	 * which case subparametric quadratic elements with six nodes are generated. The fourth, fifth, 
	 * and sixth nodes lie on the midpoints of the edges opposite the first, second, and third vertices.
	 *
	 * @param      out_tri   Output 2D triangle array. Returned by quote.
	 * @param[in]  filename  Input filename without any extensions.
	 * @param[in]  in_node   2D vertex array which contains vertex's information of the triangles.
	 * @param[in]  packed    Indicates whether the starting index of edges is zero or not. The 
	 * edge's ordering is assumed to be start with one if this option is set to NotPacked. 
	 * The default value is Packed.
	 * @param      attri     If the pointer is not null. Return a 2D integer array that contains
	 * attributes of the vertice.
	 */
	template <typename E, typename A>
	void read_Triangle_element(std::string filename, array<type_triangle2d2o<E>> &out_tri, 
		const array<vertex<point2dc, A>> &in_node, index_packed_e packed = Packed, 
		_2d_matrix *attri = nullptr)
	{
		std::ifstream tetin;
		open_infile(tetin, filename, ".ele");

		std::string err_str, tmp_str;
		std::stringstream tmp_ss;
		int ele_size = 0, node_per_ele = 0, attri_num = 0;

		// read head info
		do
		{
			std::getline(tetin, tmp_str);
			if (tmp_str[0] != '#')
			{
				str2ss(tmp_str, tmp_ss);
				tmp_ss >> ele_size >> node_per_ele >> attri_num;
				if (tmp_ss.fail() || node_per_ele != 6 || attri_num < 0)
				{
					err_str = "Wrong head-info found in " + filename + 
						".ele. From void gctl::read_Triangle_element(...)";
					throw runtime_error(err_str);
				}
			}
		}
		while (tmp_str[0] == '#');

		if (attri_num != 0 && attri != nullptr)
		{
			attri->resize(ele_size, attri_num);
		}
		else if (attri_num == 0 && attri != nullptr)
		{
			err_str = "No attributes found in " + filename + 
				".ele. From void gctl::read_Triangle_element(...)";
			throw runtime_error(err_str);
		}

		int tmp_int, tmp_order[6];
		out_tri.resize(ele_size);
		if (attri_num)
		{
			array<double> attri_line(attri_num);
			_2d_matrix tmp_attri(ele_size, attri_num);
			while (std::getline(tetin, tmp_str))
			{
				if (tmp_str[0] != '#')
				{
					str2ss(tmp_str, tmp_ss);
					tmp_ss >> tmp_int >> tmp_order[0] >> tmp_order[1] >> tmp_order[2] 
						>> tmp_order[3] >> tmp_order[4] >> tmp_order[5];
					for (int a = 0; a < attri_num; a++)
						tmp_ss >> attri_line[a];

					if (tmp_ss.fail())
					{
						err_str = "Wrong file format found in " + filename + 
							".ele. From void gctl::read_Triangle_element(...)";
						throw runtime_error(err_str);
					}

					if (packed == NotPacked)
					{
						tmp_int -= 1;
						out_tri[tmp_int].id = tmp_int;
						for (int j = 0; j < 6; j++)
						{
							out_tri[tmp_int].vert[j] = in_node.get(tmp_order[j]-1);
						}
					}
					else
					{
						out_tri[tmp_int].id = tmp_int;
						for (int j = 0; j < 6; j++)
						{
							out_tri[tmp_int].vert[j] = in_node.get(tmp_order[j]);
						}
					}

					for (int a = 0; a < attri_num; a++)
						tmp_attri[tmp_int][a] = attri_line[a];
				}
			}

			if (attri != nullptr)
			{
				for (int i = 0; i < ele_size; i++)
				{
					for (int a = 0; a < attri_num; a++)
					{
						attri->at(i, a) = tmp_attri[i][a];
					}
				}
			}
		}
		else
		{
			while (std::getline(tetin, tmp_str))
			{
				if (tmp_str[0] != '#')
				{
					str2ss(tmp_str, tmp_ss);
					tmp_ss >> tmp_int >> tmp_order[0] >> tmp_order[1] >> tmp_order[2] 
						>> tmp_order[3] >> tmp_order[4] >> tmp_order[5];

					if (tmp_ss.fail())
					{
						err_str = "Wrong file format found in " + filename + 
							".ele. From void gctl::read_Triangle_element(...)";
						throw runtime_error(err_str);
					}

					if (packed == NotPacked)
					{
						tmp_int -= 1;
						out_tri[tmp_int].id = tmp_int;
						for (int j = 0; j < 6; j++)
						{
							out_tri[tmp_int].vert[j] = in_node.get(tmp_order[j]-1);
						}
					}
					else
					{
						out_tri[tmp_int].id = tmp_int;
						for (int j = 0; j < 6; j++)
						{
							out_tri[tmp_int].vert[j] = in_node.get(tmp_order[j]);
						}
					}
				}
			}
		}

		tetin.close();
		return;
	}

	/**
	 * @brief      Read Triangle's .neigh file into an existing 2D triangle array.
	 * 
	 * First line: <# of triangles> <# of neighbors per triangle (always 3)>
	 * Following lines: <triangle #> <neighbor> <neighbor> <neighbor>
	 * Blank lines and comments prefixed by `#' may be placed anywhere. Triangles are numbered 
	 * consecutively, starting from one or zero. Neighbors are indices into the corresponding .ele file. 
	 * An index of -1 indicates no neighbor (because the triangle is on an exterior boundary). 
	 * The first neighbor of triangle i is opposite the first corner of triangle i, and so on.
	 * Triangle can produce .neigh files (use the -n switch), but cannot read them.
	 *
	 * @param      out_tri   Output 2D triangle array. Returned by quote. This array must be 
	 * initialized before calling this function.
	 * @param[in]  filename  Input filename without any extensions.
	 * @param[in]  packed    Indicates whether the starting index of edges is zero or not. The 
	 * edge's ordering is assumed to be start with one if this option is set to NotPacked. 
	 * The default value is Packed.
	 */
	template <typename A>
	void read_Triangle_neighbor(std::string filename, array<type_triangle2d<A>> &out_tri, index_packed_e packed = Packed)
	{
		std::ifstream tetin;
		open_infile(tetin, filename, ".neigh");

		std::string err_str, tmp_str;
		std::stringstream tmp_ss;
		int ele_size, neigh_pre_ele;

		// read head info
		do
		{
			std::getline(tetin, tmp_str);
			if (tmp_str[0] != '#')
			{
				str2ss(tmp_str, tmp_ss);
				tmp_ss >> ele_size >> neigh_pre_ele;
				if (tmp_ss.fail() || neigh_pre_ele != 3 || ele_size != out_tri.size())
				{
					err_str = "Wrong head-info found in " + filename + 
						".neigh. From void gctl::read_Triangle_neighbor(...)";
					throw runtime_error(err_str);
				}
			}
		}
		while (tmp_str[0] == '#');

		int tmp_int, tmp_order[3];
		while (std::getline(tetin, tmp_str))
		{
			if (tmp_str[0] != '#')
			{
				str2ss(tmp_str, tmp_ss);
				tmp_ss >> tmp_int >> tmp_order[0] >> tmp_order[1] >> tmp_order[2];

				if (tmp_ss.fail())
				{
					err_str = "Wrong file format found in " + filename + 
						".neigh. From void gctl::read_Triangle_neighbor(...)";
					throw runtime_error(err_str);
				}

				if (packed == NotPacked)
				{
					tmp_int -= 1;
					for (int j = 0; j < 3; j++)
					{
						if (tmp_order[j] > 0)
						{
							out_tri[tmp_int].neigh[j] = out_tri.get(tmp_order[j]-1);
						}
					}
				}
				else
				{
					for (int j = 0; j < 3; j++)
					{
						if (tmp_order[j] >= 0)
						{
							out_tri[tmp_int].neigh[j] = out_tri.get(tmp_order[j]);
						}
					}
				}
			}
		}

		tetin.close();
		return;
	}

	/**
	 * @brief      Read Triangle's .edge file into an 2D edge array.
	 * 
	 * First line: <# of edges> <# of boundary markers (0 or 1)>
	 * Following lines: <edge #> <endpoint> <endpoint> [boundary marker]
	 * Blank lines and comments prefixed by `#' may be placed anywhere. Edges are numbered consecutively, 
	 * starting from one or zero. Endpoints are indices into the corresponding .node file.
	 * Triangle can produce .edge files (use the -e switch), but cannot read them. The optional 
	 * column of boundary markers is suppressed by the -B switch.
	 * In Voronoi diagrams, one also finds a special kind of edge that is an infinite ray with only 
	 * one endpoint. For these edges, a different format is used:
	 * <edge #> <endpoint> -1 <direction x> <direction y>
	 * The `direction' is a floating-point vector that indicates the direction of the infinite ray.
	 *
	 * @param      out_edge  Output 2D edge array. Returned by quote.
	 * @param[in]  filename  Input filename without any extensions.
	 * @param[in]  in_node   2D vertex array which contains vertex's information of the edges.
	 * @param[in]  packed    Indicates whether the starting index of edges is zero or not. The 
	 * edge's ordering is assumed to be start with one if this option is set to NotPacked. 
	 * The default value is Packed.
	 * @param      marker   If the pointer is not null. Return a int array of the vertice's marker
	 */
	template <typename E, typename A>
	void read_Triangle_edge(std::string filename, array<type_edge2d<E>> &out_edge, 
		const array<vertex<point2dc, A>> &in_node, index_packed_e packed = Packed, 
		array<int> *marker = nullptr)
	{
		std::ifstream tetin;
		open_infile(tetin, filename, ".edge");

		std::string err_str, tmp_str;
		std::stringstream tmp_ss;
		int edge_size = 0, boundary_mark = 0;

		// read head info
		do
		{
			std::getline(tetin, tmp_str);
			if (tmp_str[0] != '#')
			{
				str2ss(tmp_str, tmp_ss);
				tmp_ss >> edge_size >> boundary_mark;
				if (tmp_ss.fail() || (boundary_mark != 0 && boundary_mark != 1))
				{
					err_str = "Wrong head-info found in " + filename + 
						".edge. From void gctl::read_Triangle_edge(...)";
					throw runtime_error(err_str);
				}
			}
		}
		while (tmp_str[0] == '#');

		if (boundary_mark != 0 && marker != nullptr)
		{
			marker->resize(edge_size);
		}
		else if (boundary_mark == 0 && marker != nullptr)
		{
			err_str = "No boundary marks found in " + filename + 
				".edge. From void gctl::read_Triangle_edge(...)";
			throw runtime_error(err_str);
		}

		int tmp_int, tmp_mark, tmp_order[2];
		out_edge.resize(edge_size);
		if (boundary_mark)
		{
			while (std::getline(tetin, tmp_str))
			{
				if (tmp_str[0] != '#')
				{
					str2ss(tmp_str, tmp_ss);
					tmp_ss >> tmp_int >> tmp_order[0] >> tmp_order[1] >> tmp_mark;

					if (tmp_ss.fail())
					{
						err_str = "Wrong file format found in " + filename + 
							".edge. From void gctl::read_Triangle_edge(...)";
						throw runtime_error(err_str);
					}

					if (packed == NotPacked)
					{
						tmp_int -= 1;
						out_edge[tmp_int].id = tmp_int;
						for (int j = 0; j < 2; j++)
						{
							out_edge[tmp_int].vert[j] = in_node.get(tmp_order[j]-1);
						}
					}
					else
					{
						out_edge[tmp_int].id = tmp_int;
						for (int j = 0; j < 2; j++)
						{
							out_edge[tmp_int].vert[j] = in_node.get(tmp_order[j]);
						}
					}

					marker->at(tmp_int) = tmp_mark;
				}
			}
		}
		else
		{
			while (std::getline(tetin, tmp_str))
			{
				if (tmp_str[0] != '#')
				{
					str2ss(tmp_str, tmp_ss);
					tmp_ss >> tmp_int >> tmp_order[0] >> tmp_order[1];

					if (tmp_ss.fail())
					{
						err_str = "Wrong file format found in " + filename + 
							".edge. From void gctl::read_Triangle_edge(...)";
						throw runtime_error(err_str);
					}

					if (packed == NotPacked)
					{
						tmp_int -= 1;
						out_edge[tmp_int].id = tmp_int;
						for (int j = 0; j < 2; j++)
						{
							out_edge[tmp_int].vert[j] = in_node.get(tmp_order[j]-1);
						}
					}
					else
					{
						out_edge[tmp_int].id = tmp_int;
						for (int j = 0; j < 2; j++)
						{
							out_edge[tmp_int].vert[j] = in_node.get(tmp_order[j]);
						}
					}
				}
			}
		}

		tetin.close();
		return;
	}

	/**
	 * @brief      Saves a Triangle node file.
	 *
	 * @param[in]  filename  Save name
	 * @param[in]  out_node  Output 2D vertice
	 * @param[in]  packed    Indicates whether the starting index of edges is zero or not. The 
	 * edge's ordering is assumed to be start with one if this option is set to NotPacked. 
	 * The default value is Packed.
	 * @param[in]  marker    If the pointer is not null. Write a int array of the vertice's marker
	 * @param      attri     If the pointer is not null. Write a 2D integer array that contains
	 * attributes of the vertice.
	 */
	template <typename A>
	void save_Triangle_node(std::string filename, const array<vertex<point2dc, A>> &out_node, 
		index_packed_e packed = Packed, const array<int> *marker = nullptr, const _2d_matrix *attri = nullptr)
	{
		std::ofstream triout;
		open_outfile(triout, filename, ".node");

		std::string err_str;
		int attri_num = 0, boundary_mark = 0;
		if (marker != nullptr)
		{
			boundary_mark = 1;
			if (out_node.size() != marker->size())
			{
				err_str = "Size of the boundary marks does not match. From void gctl::save_Triangle_node(...)";
				throw runtime_error(err_str);
			}
		}

		if (attri != nullptr)
		{
			attri_num = attri->col_size();
			if (out_node.size() != attri->row_size())
			{
				err_str = "Size of the attributes do not match. From void gctl::save_Triangle_node(...)";
				throw runtime_error(err_str);
			}
		}

		time_t now = time(0);
		char *dt = ctime(&now);
		triout << "# generated by the GCTL package on " << dt;
		triout << "# node_num node_dimension attri_num boundary_mark" << std::endl;
		triout << out_node.size() << " 2 " << attri_num << " " << boundary_mark << std::endl;
		triout << "# node_index x y attributes mark" << std::endl;
		if (attri != nullptr && marker != nullptr)
		{
			if (packed == Packed)
			{
				for (int i = 0; i < out_node.size(); i++)
				{
					triout << i << " " << std::setprecision(16) << out_node[i].x << " " 
						<< out_node[i].y;

					for (int j = 0; j < attri_num; j++)
					{
						triout << " " << attri->at(i, j);
					}

					triout << " " << marker->at(i) << std::endl;
				}
			}
			else
			{
				for (int i = 0; i < out_node.size(); i++)
				{
					triout << i+1 << " " << std::setprecision(16) << out_node[i].x << " " 
						<< out_node[i].y;

					for (int j = 0; j < attri_num; j++)
					{
						triout << " " << attri->at(i, j);
					}

					triout << " " << marker->at(i) << std::endl;
				}
			}
		}
		else if (attri != nullptr)
		{
			if (packed == Packed)
			{
				for (int i = 0; i < out_node.size(); i++)
				{
					triout << i << " " << std::setprecision(16) << out_node[i].x << " " 
						<< out_node[i].y;

					for (int j = 0; j < attri_num; j++)
					{
						triout << " " << attri->at(i, j);
					}

					triout << std::endl;
				}
			}
			else
			{
				for (int i = 0; i < out_node.size(); i++)
				{
					triout << i+1 << " " << std::setprecision(16) << out_node[i].x << " " 
						<< out_node[i].y;

					for (int j = 0; j < attri_num; j++)
					{
						triout << " " << attri->at(i, j);
					}

					triout << std::endl;
				}
			}
		}
		else if (marker != nullptr)
		{
			if (packed == Packed)
			{
				for (int i = 0; i < out_node.size(); i++)
				{
					triout << i << " " << std::setprecision(16) << out_node[i].x << " " 
						<< out_node[i].y << " " << marker->at(i) << std::endl;
				}
			}
			else
			{
				for (int i = 0; i < out_node.size(); i++)
				{
					triout << i+1 << " " << std::setprecision(16) << out_node[i].x << " " 
						<< out_node[i].y << " " << marker->at(i) << std::endl;
				}
			}
		}
		else
		{
			if (packed == Packed)
			{
				for (int i = 0; i < out_node.size(); i++)
				{
					triout << i << " " << std::setprecision(16) << out_node[i].x << " " 
						<< out_node[i].y << std::endl;
				}
			}
			else
			{
				for (int i = 0; i < out_node.size(); i++)
				{
					triout << i+1 << " " << std::setprecision(16) << out_node[i].x << " " 
						<< out_node[i].y << std::endl;
				}
			}
		}

		triout.close();
		return;
	}

	/**
	 * @brief      Saves a Triangle element file.
	 *
	 * @param[in]  filename  Save name.
	 * @param[in]  out_tri   Output 2D triangles.
	 * @param[in]  packed    Indicates whether the starting index of edges is zero or not. The 
	 * edge's ordering is assumed to be start with one if this option is set to NotPacked. 
	 * The default value is Packed.
	 * @param      attri     If the pointer is not null. Write a 2D integer array that contains
	 * attributes of the vertice.
	 */
	template <typename A>
	void save_Triangle_element(std::string filename, const array<type_triangle2d<A>> &out_tri, 
		index_packed_e packed = Packed, const _2d_matrix *attri = nullptr)
	{
		std::ofstream triout;
		open_outfile(triout, filename, ".ele");

		std::string err_str;
		int attri_num = 0;
		if (attri != nullptr)
		{
			attri_num = attri->col_size();
			if (out_tri.size() != attri->row_size())
			{
				err_str = "Size of the attributes do not match. From void gctl::save_Triangle_element(...)";
				throw runtime_error(err_str);
			}
		}

		time_t now = time(0);
		char *dt = ctime(&now);
		triout << "# generated by the GCTL package on " << dt;
		triout << "# element_num index_num attri_num" << std::endl;
		triout << out_tri.size() << " 3 " << attri_num << std::endl;
		triout << "# element_index node_index attributes" << std::endl;
		if (attri != nullptr)
		{
			if (packed == Packed)
			{
				for (int i = 0; i < out_tri.size(); i++)
				{
					triout << i;
					for (int j = 0; j < 3; j++)
					{
						triout << " " << out_tri[i].vert[j]->id;
					}

					for (int j = 0; j < attri_num; j++)
					{
						triout << " " << attri->at(i, j);
					}
					triout << std::endl;
				}
			}
			else
			{
				for (int i = 0; i < out_tri.size(); i++)
				{
					triout << i + 1;
					for (int j = 0; j < 3; j++)
					{
						triout << " " << out_tri[i].vert[j]->id + 1;
					}

					for (int j = 0; j < attri_num; j++)
					{
						triout << " " << attri->at(i, j);
					}
					triout << std::endl;
				}
			}
		}
		else
		{
			if (packed == Packed)
			{
				for (int i = 0; i < out_tri.size(); i++)
				{
					triout << i;
					for (int j = 0; j < 3; j++)
					{
						triout << " " << out_tri[i].vert[j]->id;
					}
					triout << std::endl;
				}
			}
			else
			{
				for (int i = 0; i < out_tri.size(); i++)
				{
					triout << i + 1;
					for (int j = 0; j < 3; j++)
					{
						triout << " " << out_tri[i].vert[j]->id + 1;
					}
					triout << std::endl;
				}
			}
		}

		triout.close();
		return;
	}

	/**
	 * @brief      Saves a Triangle element file.
	 *
	 * @param[in]  filename  Save name.
	 * @param[in]  out_tri   Output 2D triangles.
	 * @param[in]  packed    Indicates whether the starting index of edges is zero or not. The 
	 * edge's ordering is assumed to be start with one if this option is set to NotPacked. 
	 * The default value is Packed.
	 * @param      attri     If the pointer is not null. Write a 2D integer array that contains
	 * attributes of the vertice.
	 */
	template <typename A>
	void save_Triangle_element(std::string filename, const array<type_triangle2d2o<A>> &out_tri, 
		index_packed_e packed = Packed, const _2d_matrix *attri = nullptr)
	{
		std::ofstream triout;
		open_outfile(triout, filename, ".ele");

		std::string err_str;
		int attri_num = 0;
		if (attri != nullptr)
		{
			attri_num = attri->col_size();
			if (out_tri.size() != attri->row_size())
			{
				err_str = "Size of the attributes do not match. From void gctl::save_Triangle_element(...)";
				throw runtime_error(err_str);
			}
		}

		time_t now = time(0);
		char *dt = ctime(&now);
		triout << "# generated by the GCTL package on " << dt;
		triout << "# element_num index_num attri_num" << std::endl;
		triout << out_tri.size() << " 6 " << attri_num << std::endl;
		triout << "# element_index node_index attributes" << std::endl;
		if (attri != nullptr)
		{
			if (packed == Packed)
			{
				for (int i = 0; i < out_tri.size(); i++)
				{
					triout << i;
					for (int j = 0; j < 6; j++)
					{
						triout << " " << out_tri[i].vert[j]->id;
					}

					for (int j = 0; j < attri_num; j++)
					{
						triout << " " << attri->at(i, j);
					}
					triout << std::endl;
				}
			}
			else
			{
				for (int i = 0; i < out_tri.size(); i++)
				{
					triout << i + 1;
					for (int j = 0; j < 6; j++)
					{
						triout << " " << out_tri[i].vert[j]->id + 1;
					}

					for (int j = 0; j < attri_num; j++)
					{
						triout << " " << attri->at(i, j);
					}
					triout << std::endl;
				}
			}
		}
		else
		{
			if (packed == Packed)
			{
				for (int i = 0; i < out_tri.size(); i++)
				{
					triout << i;
					for (int j = 0; j < 6; j++)
					{
						triout << " " << out_tri[i].vert[j]->id;
					}
					triout << std::endl;
				}
			}
			else
			{
				for (int i = 0; i < out_tri.size(); i++)
				{
					triout << i + 1;
					for (int j = 0; j < 6; j++)
					{
						triout << " " << out_tri[i].vert[j]->id + 1;
					}
					triout << std::endl;
				}
			}
		}

		triout.close();
		return;
	}

	/**
	 * @brief      Saves a Triangle neighbor file.
	 *
	 * @param[in]  filename  Save name.
	 * @param[in]  out_tri   Output 2D triangles.
	 * @param[in]  packed    Indicates whether the starting index of edges is zero or not. The 
	 * edge's ordering is assumed to be start with one if this option is set to NotPacked. 
	 * The default value is Packed.
	 */
	template <typename A>
	void save_Triangle_neighbor(std::string filename, const array<type_triangle2d<A>> &out_tri, 
		index_packed_e packed = Packed)
	{
		std::ofstream triout;
		open_outfile(triout, filename, ".neigh");

		time_t now = time(0);
		char *dt = ctime(&now);
		triout << "# generated by the GCTL package on " << dt;
		triout << "# element_num index_num" << std::endl;
		triout << out_tri.size() << " 3" << std::endl;
		triout << "# element_index element_index" << std::endl;
		if (packed == Packed)
		{
			for (int i = 0; i < out_tri.size(); i++)
			{
				triout << i;
				for (int j = 0; j < 3; j++)
				{
					if (out_tri[i].neigh[j] != nullptr)
						triout << " " << out_tri[i].neigh[j]->id;
					else triout << " -1";
				}
				triout << std::endl;
			}
		}
		else
		{
			for (int i = 0; i < out_tri.size(); i++)
			{
				triout << i + 1;
				for (int j = 0; j < 3; j++)
				{
					if (out_tri[i].neigh[j] != nullptr)
						triout << " " << out_tri[i].neigh[j]->id + 1;
					else triout << " -1";
				}
				triout << std::endl;
			}
		}

		triout.close();
		return;
	}

	/**
	 * @brief      Saves a Triangle edge file.
	 *
	 * @param[in]  filename  Save name.
	 * @param[in]  out_edge  Output 2D edges.
	 * @param[in]  packed    Indicates whether the starting index of edges is zero or not. The 
	 * edge's ordering is assumed to be start with one if this option is set to NotPacked. 
	 * The default value is Packed.
	 * @param      bound_mark   If the pointer is not null. Write a int array of the vertice's marker
	 */
	template <typename A>
	void save_Triangle_edge(std::string filename, const array<type_edge2d<A>> &out_edge, 
		index_packed_e packed = Packed, const array<int> *marker = nullptr)
	{
		std::ofstream triout;
		open_outfile(triout, filename, ".edge");

		std::string err_str;
		int boundary_mark = 0;
		if (marker != nullptr)
		{
			boundary_mark = 1;
			if (out_edge.size() != marker->size())
			{
				err_str = "Size of the boundary marks does not match. From void gctl::save_Triangle_edge(...)";
				throw runtime_error(err_str);
			}
		}

		time_t now = time(0);
		char *dt = ctime(&now);
		triout << "# generated by the GCTL package on " << dt;
		triout << "# edge_num index_num" << std::endl;
		triout << out_edge.size() << " " << boundary_mark << std::endl;
		triout << "# edge_index node_index boundary_mark" << std::endl;
		if (marker != nullptr)
		{
			if (packed == Packed)
			{
				for (int i = 0; i < out_edge.size(); i++)
				{
					triout << i;
					for (int j = 0; j < 2; j++)
					{
						triout << " " << out_edge[i].vert[j]->id;
					}

					triout << " " << marker->at(i) << std::endl;
				}
			}
			else
			{
				for (int i = 0; i < out_edge.size(); i++)
				{
					triout << i + 1;
					for (int j = 0; j < 2; j++)
					{
						triout << " " << out_edge[i].vert[j]->id + 1;
					}

					triout << " " << marker->at(i) << std::endl;
				}
			}
		}
		else
		{
			if (packed == Packed)
			{
				for (int i = 0; i < out_edge.size(); i++)
				{
					triout << i;
					for (int j = 0; j < 2; j++)
					{
						triout << " " << out_edge[i].vert[j]->id;
					}
					triout << std::endl;
				}
			}
			else
			{
				for (int i = 0; i < out_edge.size(); i++)
				{
					triout << i + 1;
					for (int j = 0; j < 2; j++)
					{
						triout << " " << out_edge[i].vert[j]->id + 1;
					}
					triout << std::endl;
				}
			}
		}

		triout.close();
		return;
	}

};

#endif //_GCTL_TRIANGLE_H