/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * 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_TETGEN_IO_H
#define _GCTL_TETGEN_IO_H

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

namespace gctl
{
	/**
	 * @brief      Read Tetrahedron's .node file into a 3D vertex array.
	 * 
	 * First line:  <# of points> <dimension (3)> <# of attributes> <boundary markers (0 or 1)>
	 * Remaining lines list # of points: 
	 * <point #> <x> <y> <z> [attributes] [boundary marker]
	 * All tetgen files are of ASCII form and may contain comments prefixed by the 
	 * character ’#’. Points, tetrahedra, facets, edges, holes, and maximum volume constraints must 
	 * be numbered consecutively, starting from either 1 or 0. However, all input files must be 
	 * consistent. TetGen automatically detects your choice while reading the .node 
	 * (or .poly or .smesh) file. When calling TetGen from another program, 
	 * use the -z switch if you wish to number objects from zero.
	 * 
	 * @param[in]  out_node  Output 3D vertex array. Returned by quote.
	 * @param      filename  Input filename without any extensions.
	 * @param[in]  packed    Indicates whether the starting index of vertice is zero or not. The 
	 * vertex's ordering is assumed to be start with one if this option is set to false. 
	 * The default value is true.
	 * @param[in]  bound_tag Pointer of the boundary marker array.
	 * @param[in]  attri     Pointer of the attributes array.
	 */
	template <typename A>
	void read_Tetgen_node(std::string filename, array<vertex<point3dc, A>> &out_node, index_packed_e packed = Packed, 
		array<int> *bound_tag = 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 = 3, 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 != 3 || 
					(boundary_mark != 0 && boundary_mark != 1) || 
					attri_num < 0)
				{
					err_str = "Wrong head-info found in " + filename + 
						".node. From void gctl::read_Tetgen_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_Tetgen_node(...)";
			throw runtime_error(err_str);
		}

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

		vertex<point3dc, 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 >> tmp_vert.z;
					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_Tetgen_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)
			{
				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 (bound_tag != nullptr)
			{
				for (int i = 0; i < node_size; i++)
				{
					bound_tag->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 >> tmp_vert.z;

					if (tmp_ss.fail())
					{
						err_str = "Wrong file format found in " + filename + 
							".node. From void gctl::read_Tetgen_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 Tetgen's .ele file into a 3D tetrahedron array.
	 * 
	 * First line: <# of tetrahedra> <nodes per tet. (4 or 10)> <region attribute (0 or 1)>
	 * Remaining lines list # of tetrahedra:
	 * <tetrahedron #> <node> <node> ... <node> [attribute]
	 *
	 * @param      out_tri   Output 3D tetrahedron array. Returned by quote.
	 * @param[in]  filename  Input filename without any extensions.
	 * @param[in]  in_node   3D vertex array which contains vertex's information of the tetrahedrons.
	 * @param[in]  packed    Indicates whether the starting index of vertice is zero or not. The 
	 * vertex's ordering is assumed to be start with one if this option is set to false. 
	 * The default value is true.
	 * @param[in]  region_attri  If the <region attribute> in the first line is 1, 
	 * each tetrahedra has an additional region attribute (an integer) in the last column. 
	 */
	template <typename E, typename A>
	void read_Tetgen_element(std::string filename, array<type_tetrahedron<E>> &out_tet, 
		const array<vertex<point3dc, A>> &in_node, index_packed_e packed = Packed, 
		array<int> *region_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 = 4, region_mark = 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 >> region_mark;
				if (tmp_ss.fail() || node_per_ele != 4 || 
					(region_mark != 0 && region_mark != 1))
				{
					err_str = "Wrong head-info found in " + filename + 
						".ele. From void gctl::read_Tetgen_element(...)";
					throw runtime_error(err_str);
				}
			}
		}
		while (tmp_str[0] == '#');

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

		int tmp_int, tmp_order[4], tmp_region;
		out_tet.resize(ele_size);
		if (region_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_order[2] 
						>> tmp_order[3] >> tmp_region;

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

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

						if (region_attri != nullptr)
						{
							region_attri->at(tmp_int) = tmp_region;
						}
					}
					else
					{
						out_tet[tmp_int].id = tmp_int;
						for (int j = 0; j < 4; j++)
						{
							out_tet[tmp_int].vert[j] = in_node.get(tmp_order[j]);
						}
						out_tet[tmp_int].deter_vert_order();

						if (region_attri != nullptr)
						{
							region_attri->at(tmp_int) = tmp_region;
						}
					}
				}
			}
		}
		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];

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

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

		tetin.close();
		return;
	}

	/**
	 * @brief       Read neighbor file generated by the Tetgen program.
	 * 
	 * First line:  <# of tetrahedra>  4
	 * Following lines list # of neighbors:
	 * <tetrahedra #> <neighbor> <neighbor> <neighbor> <neighbor>
	 * An index of −1 indicates no neighbor (because the tetrahedron is 
	 * on boundary of a mesh domain).
	 *
	 * @param      out_tet   The output tetrahedron array.
	 * @param[in]  filename  The filename of a .neigh file generated by the Tetgen program.
	 * @param[in]  packed    Indicates whether the index in the neighbor file starts from zero. The 
	 * index is deemed to be started with one if this option is false. The default value of this 
	 * variable is true.
	 */
	template <typename A>
	void read_Tetgen_neighbor(std::string filename, array<type_tetrahedron<A>> &out_tet, 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 != 4 || ele_size != out_tet.size())
				{
					err_str = "Wrong head-info found in " + filename + 
						".neigh. From void gctl::read_Tetgen_neighbor(...)";
					throw runtime_error(err_str);
				}
			}
		}
		while (tmp_str[0] == '#');

		int tmp_int, tmp_order[4];
		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];

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

				if (packed == NotPacked)
				{
					tmp_int -= 1;
					for (int j = 0; j < 4; j++)
					{
						if (tmp_order[j] != -1) out_tet[tmp_int].neigh[j] = out_tet.get(tmp_order[j]-1);
						else out_tet[tmp_int].neigh[j] = nullptr;
					}
				}
				else
				{
					for (int j = 0; j < 4; j++)
					{
						if (tmp_order[j] != -1) out_tet[tmp_int].neigh[j] = out_tet.get(tmp_order[j]);
						else out_tet[tmp_int].neigh[j] = nullptr;
					}
				}
			}
		}

		tetin.close();
		return;
	}

	/**
	 * @brief      Read edge file generated by the Tetgen program.
	 * 
	 * First line: <# of edges> <boundary marker (0 or 1)>
	 * Remaining lines list # of edges:
	 * <edge #> <endpoint> <endpoint> [boundary marker]
	 *
	 * @param      out_edge  The output edge array.
	 * @param[in]  filename  The filename of a .edge file generated by the Tetgen program.
	 * @param[in]  in_node   3D vertex array which contains vertex's information of the tetrahedrons.
	 * @param[in]  packed    Indicates whether the starting index of vertice is zero or not. The 
	 * vertex's ordering is assumed to be start with one if this option is set to false. 
	 * The default value is true.
	 * @param[in]  bound_tag Pointer of the boundary marker array.
	 */
	template <typename E, typename A>
	void read_Tetgen_edge(std::string filename, array<type_edge<E>> &out_edge, 
		const array<vertex<point3dc, A>> &in_node, index_packed_e packed = Packed, 
		array<int> *bound_tag = 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_Tetgen_edge(...)";
					throw runtime_error(err_str);
				}
			}
		}
		while (tmp_str[0] == '#');

		if (boundary_mark != 0 && bound_tag != nullptr)
		{
			bound_tag->resize(edge_size);
		}
		else if (boundary_mark == 0 && bound_tag != nullptr)
		{
			err_str = "No boundary marks found in " + filename + 
				".edge. From void gctl::read_Tetgen_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_Tetgen_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]);
						}
					}

					bound_tag->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_Tetgen_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       Read face file generated by the Tetgen program.
	 * 
	 * First line: <# of faces> <boundary marker (0 or 1)>
	 * Remaining lines list # of faces:
	 * <face #> <node> <node> <node> ... [boundary marker] ...
	 *
	 * @param      out_tri    The output face array.
	 * @param[in]  filename   The filename of a .face file generated by the Tetgen program.
	 * @param[in]  in_node    3D vertex array which contains vertex's information of the tetrahedrons.
	 * @param[in]  packed     Indicates whether the starting index of vertice is zero or not. The 
	 * vertex's ordering is assumed to be start with one if this option is set to false. 
	 * The default value is true.
	 * @param[in]  bound_tag Pointer of the boundary marker array.
	 */
	template <typename E, typename A>
	void read_Tetgen_face(std::string filename, array<type_triangle<E>> &out_tri, 
		const array<vertex<point3dc, A>> &in_node, index_packed_e packed = Packed, 
		array<int> *bound_tag = nullptr)
	{
		std::ifstream tetin;
		open_infile(tetin, filename, ".face");

		std::string err_str, tmp_str;
		std::stringstream tmp_ss;
		int face_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 >> face_size >> boundary_mark;
				if (tmp_ss.fail() || (boundary_mark != 0 && boundary_mark != 1))
				{
					err_str = "Wrong head-info found in " + filename + 
						".face. From void gctl::read_Tetgen_face(...)";
					throw runtime_error(err_str);
				}
			}
		}
		while (tmp_str[0] == '#');

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

		int tmp_int, tmp_mark, tmp_order[3];
		out_tri.resize(face_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_order[2] >> tmp_mark;

					if (tmp_ss.fail())
					{
						err_str = "Wrong file format found in " + filename + 
							".face. From void gctl::read_Tetgen_face(...)";
						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]);
						}
					}

					if (bound_tag != nullptr)
					{
						bound_tag->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] 
						>> tmp_order[2];

					if (tmp_ss.fail())
					{
						err_str = "Wrong file format found in " + filename + 
							".face. From void gctl::read_Tetgen_face(...)";
						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      Saves a Tetgen node file.
	 *
	 * @param[in]  filename  Output filename.
	 * @param[in]  out_node  Output nodes.
	 * @param[in]  packed    Indicates whether the starting index of vertice is zero or not. The 
	 * vertex's ordering is assumed to be start with one if this option is set to false. 
	 * The default value is true.
	 */
	template <typename A>
	void save_Tetgen_node(std::string filename, const array<vertex<point3dc, A>> &out_node, index_packed_e packed = Packed, 
		const array<int> *bound_tag = 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 (bound_tag != nullptr)
		{
			boundary_mark = 1;
			if (out_node.size() != bound_tag->size())
			{
				err_str = "Size of the boundary marks does not match. From void gctl::save_Tetgen_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_Tetgen_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() << " 3 " << attri_num << " " << boundary_mark << std::endl;
		triout << "# node_index x y z attributes mark" << std::endl;
		if (attri != nullptr && bound_tag != 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 << " " << out_node[i].z;

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

					triout << " " << bound_tag->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 << " " << out_node[i].z;

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

					triout << " " << bound_tag->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 << " " << out_node[i].z;

					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 << " " << out_node[i].z;

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

					triout << std::endl;
				}
			}
		}
		else if (bound_tag != 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 << " " << out_node[i].z;

					triout << " " << bound_tag->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 << " " << out_node[i].z;

					triout << " " << bound_tag->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 << " " << out_node[i].z << 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 << " " << out_node[i].z << std::endl;
				}
			}
		}

		triout.close();
		return;
	}

	/**
	 * @brief      Saves a Tetgen element file.
	 *
	 * @param[in]  filename  Output filename.
	 * @param[in]  out_tet   Output tetrahedral elements.
	 * @param[in]  packed    Indicates whether the starting index of vertice is zero or not. The 
	 * vertex's ordering is assumed to be start with one if this option is set to false. 
	 * The default value is true.
	 */
	template <typename A>
	void save_Tetgen_element(std::string filename, const array<type_tetrahedron<A>> &out_tet, 
		index_packed_e packed = Packed, const array<int> *region_attri = nullptr)
	{
		std::ofstream tetout;
		open_outfile(tetout, filename, ".ele");

		std::string err_str;
		int region_mark = 0;
		if (region_attri != nullptr)
		{
			region_mark = 1;
			if (out_tet.size() != region_attri->size())
			{
				err_str = "Size of the region attribute does not match. From void gctl::save_Tetgen_element(...)";
				throw runtime_error(err_str);
			}
		}

		time_t now = time(0);
		char *dt = ctime(&now);
		tetout << "# generated by the GCTL package on " << dt;
		tetout << "# element_num index_num region_mark" << std::endl;
		tetout << out_tet.size() << " 4 " << region_mark << std::endl;
		tetout << "# element_index node_index attributes" << std::endl;
		if (region_attri != nullptr)
		{
			if (packed == Packed)
			{
				for (int i = 0; i < out_tet.size(); i++)
				{
					tetout << i;
					for (int j = 0; j < 4; j++)
					{
						tetout << " " << out_tet[i].vert[j]->id;
					}
					tetout << " " << region_attri->at(i) << std::endl;
				}
			}
			else
			{
				for (int i = 0; i < out_tet.size(); i++)
				{
					tetout << i + 1;
					for (int j = 0; j < 4; j++)
					{
						tetout << " " << out_tet[i].vert[j]->id + 1;
					}
					tetout << " " << region_attri->at(i) << std::endl;
				}
			}
		}
		else
		{
			if (packed == Packed)
			{
				for (int i = 0; i < out_tet.size(); i++)
				{
					tetout << i;
					for (int j = 0; j < 4; j++)
					{
						tetout << " " << out_tet[i].vert[j]->id;
					}
					tetout << std::endl;
				}
			}
			else
			{
				for (int i = 0; i < out_tet.size(); i++)
				{
					tetout << i + 1;
					for (int j = 0; j < 4; j++)
					{
						tetout << " " << out_tet[i].vert[j]->id + 1;
					}
					tetout << std::endl;
				}
			}
		}

		tetout.close();
		return;
	}

	/**
	 * @brief      Saves a Tetgen neighbor file.
	 *
	 * @param[in]  filename  Output filename.
	 * @param[in]  out_tet   Output tetrahedral elements.
	 * @param[in]  packed    Indicates whether the starting index of vertice is zero or not. The 
	 * vertex's ordering is assumed to be start with one if this option is set to false. 
	 * The default value is true.
	 */
	template <typename A>
	void save_Tetgen_neighbor(std::string filename, const array<type_tetrahedron<A>> &out_tet, 
		index_packed_e packed = Packed)
	{
		std::ofstream tetout;
		open_outfile(tetout, filename, ".neigh");

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

		tetout.close();
		return;
	}

	/**
	 * @brief      Saves a Tetgen edge file.
	 *
	 * @param[in]  filename  Output filename.
	 * @param      out_edge  Output edge elements.
	 * @param[in]  packed    Indicates whether the starting index of vertice is zero or not. The 
	 * vertex's ordering is assumed to be start with one if this option is set to false. 
	 * The default value is true.
	 */
	template <typename A>
	void save_Tetgen_edge(std::string filename, const array<type_edge<A>> &out_edge, index_packed_e packed = Packed, 
		const array<int> *bound_tag = nullptr)
	{
		std::ofstream tetout;
		open_outfile(tetout, filename, ".edge");

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

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

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

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

		tetout.close();
		return;
	}

	/**
	 * @brief      Saves a Tetgen face face.
	 *
	 * @param[in]  filename  Output filename.
	 * @param[in]  out_tri   Output triangle elements.
	 * @param[in]  packed    Indicates whether the starting index of vertice is zero or not. The 
	 * vertex's ordering is assumed to be start with one if this option is set to false. 
	 * The default value is true.
	 */
	template <typename A>
	void save_Tetgen_face(std::string filename, const array<type_triangle<A>> &out_tri, index_packed_e packed = Packed, 
		const array<int> *bound_tag = nullptr)
	{
		std::ofstream tetout;
		open_outfile(tetout, filename, ".face");

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

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

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

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

		tetout.close();
		return;
	}
};

#endif //_GCTL_TETGEN_IO_H