/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * 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.
 ******************************************************/

#include "regular_grid.h"

#ifdef GCTL_MESH_WAVELIB
// include wavelib headfile
#include "wavelib.h"
#endif // GCTL_MESH_WAVELIB

void gctl::regular_grid::init(std::string in_name, std::string in_info, int xnum, int ynum, 
	double xmin, double ymin, double dx, double dy)
{
	if (initialized)
	{
		throw std::runtime_error("[gctl::regular_grid] The mesh is already initialized.");
	}

	base_mesh::init(REGULAR_GRID, MESH_2D, in_name, in_info);

	rg_xnum = xnum; rg_ynum = ynum;
	rg_xmin = xmin; rg_ymin = ymin;
	rg_dx = dx; rg_dy = dy;

	node_num = rg_xnum * rg_ynum;
	ele_num  = (rg_xnum-1) * (rg_ynum-1);

	nodes.resize(node_num);
	elements.resize(ele_num);

	for (int j = 0; j < rg_ynum; j++)
	{
		for (int i = 0; i < rg_xnum; i++)
		{
			nodes[i+j*rg_xnum].id = i + j*rg_xnum;
			nodes[i+j*rg_xnum].x  = rg_xmin + i*rg_dx;
			nodes[i+j*rg_xnum].y  = rg_ymin + j*rg_dy;
		}
	}

	for (int j = 0; j < rg_ynum-1; j++)
	{
		for (int i = 0; i < rg_xnum-1; i++)
		{
			elements[i+j*(rg_xnum-1)].id = i+j*(rg_xnum-1);
			elements[i+j*(rg_xnum-1)].vert[0] = nodes.get(i+j*rg_xnum);
			elements[i+j*(rg_xnum-1)].vert[1] = nodes.get(i+1+j*rg_xnum);
			elements[i+j*(rg_xnum-1)].vert[2] = nodes.get(i+1+(j+1)*rg_xnum);
			elements[i+j*(rg_xnum-1)].vert[3] = nodes.get(i+(j+1)*rg_xnum);
			elements[i+j*(rg_xnum-1)].dl = nodes.get(i+j*rg_xnum);
			elements[i+j*(rg_xnum-1)].ur = nodes.get(i+1+(j+1)*rg_xnum);
		}
	}

	initialized = true;
	return;
}

void gctl::regular_grid::show_mesh_dimension(std::ostream &os) const
{
	os << "Range X: " << rg_xmin << " -> " << rg_xmin + (rg_xnum - 1)*rg_dx << "\n";
	os << "Range Y: " << rg_ymin << " -> " << rg_ymin + (rg_ynum - 1)*rg_dy << "\n";
	os << "Interval: " << rg_dx << ", " << rg_dy << "\n";
	return;
}

void gctl::regular_grid::load_binary(std::string filename)
{
	if (initialized)
	{
		throw std::runtime_error("[gctl::regular_grid] The mesh is already initialized.");
	}

	std::ifstream infile;
	gctl::open_infile(infile, filename, ".2m", std::ios::in|std::ios::binary);

	// 读入网格头信息
	load_headinfo(infile, REGULAR_GRID, MESH_2D);

	// 读入网格信息
	infile.read((char*)&rg_xnum, sizeof(int));
	infile.read((char*)&rg_ynum, sizeof(int));
	infile.read((char*)&rg_xmin, sizeof(double));
	infile.read((char*)&rg_ymin, sizeof(double));
	infile.read((char*)&rg_dx, sizeof(double));
	infile.read((char*)&rg_dy, sizeof(double));

	node_num = rg_xnum * rg_ynum;
	ele_num  = (rg_xnum-1) * (rg_ynum-1);

	nodes.resize(node_num);
	elements.resize(ele_num);

	for (int j = 0; j < rg_ynum; j++)
	{
		for (int i = 0; i < rg_xnum; i++)
		{
			nodes[i+j*rg_xnum].id = i + j*rg_xnum;
			nodes[i+j*rg_xnum].x  = rg_xmin + i*rg_dx;
			nodes[i+j*rg_xnum].y  = rg_ymin + j*rg_dy;
		}
	}

	for (int j = 0; j < rg_ynum-1; j++)
	{
		for (int i = 0; i < rg_xnum-1; i++)
		{
			elements[i+j*(rg_xnum-1)].id = i+j*(rg_xnum-1);
			elements[i+j*(rg_xnum-1)].vert[0] = nodes.get(i+j*rg_xnum);
			elements[i+j*(rg_xnum-1)].vert[1] = nodes.get(i+1+j*rg_xnum);
			elements[i+j*(rg_xnum-1)].vert[2] = nodes.get(i+1+(j+1)*rg_xnum);
			elements[i+j*(rg_xnum-1)].vert[3] = nodes.get(i+(j+1)*rg_xnum);
			elements[i+j*(rg_xnum-1)].dl = nodes.get(i+j*rg_xnum);
			elements[i+j*(rg_xnum-1)].ur = nodes.get(i+1+(j+1)*rg_xnum);
		}
	}

	initialized = true;

	// 读入模型数据单元
	load_datablock(infile);
	infile.close();	
	return;
}

void gctl::regular_grid::save_binary(std::string filename)
{
	if (!initialized)
	{
		throw std::runtime_error("[gctl::regular_grid] The mesh is not initialized.");
	}

	std::ofstream outfile;
	gctl::open_outfile(outfile, filename, ".2m", std::ios::out|std::ios::binary);

	// 首先输出网格的头信息
	save_headinfo(outfile);

	// 输出网格信息
	outfile.write((char*)&rg_xnum, sizeof(int));
	outfile.write((char*)&rg_ynum, sizeof(int));
	outfile.write((char*)&rg_xmin, sizeof(double));
	outfile.write((char*)&rg_ymin, sizeof(double));
	outfile.write((char*)&rg_dx, sizeof(double));
	outfile.write((char*)&rg_dy, sizeof(double));

	// 输出的模型数据单元
	save_datablock(outfile);

	outfile.close();
	return;
}

gctl::regular_grid::regular_grid() : base_mesh::base_mesh(){}

gctl::regular_grid::regular_grid(std::string in_name, std::string in_info, int xnum, int ynum, 
	double xmin, double ymin, double dx, double dy)
{
	init(in_name, in_info, xnum, ynum, xmin, ymin, dx, dy);
}

gctl::regular_grid::~regular_grid(){}

void gctl::regular_grid::clear_regular_grid()
{
	clear();
	elements.clear();
	nodes.clear();
	return;
}

#ifdef GCTL_GRAPHIC_MATHGL

int gctl::regular_grid::view(std::string datname)
{
	array<double> *datval_ptr = (array<double>*) get_datval(datname);

	plt.range(rg_xmin, rg_xmin + (rg_xnum - 1)*rg_dx, rg_ymin, rg_ymin + (rg_ynum - 1)*rg_dy);
	plt.demension(rg_xnum, rg_ynum);
	plt.add_dens(*datval_ptr, datname);
	return plt.plot();
}

#else

int gctl::regular_grid::view(std::string datname)
{
	return 0;
}

#endif // GCTL_GRAPHIC_MATHGL

#ifdef GCTL_GRAPHIC_GMT

void gctl::regular_grid::plot(std::string datname)
{
	array<double> *datval_ptr = (array<double>*) get_datval(datname);

	pic.set_command("psconvert", "-A -TG -E300");
	pic.plot(datname, *datval_ptr, 
		rg_xmin, rg_xmin + (rg_xnum - 1)*rg_dx, 
		rg_ymin, rg_ymin + (rg_ynum - 1)*rg_dy,
		rg_xnum, rg_ynum);
	return;
}

#else

void gctl::regular_grid::plot(std::string datname)
{
	return;
}

#endif // GCTL_GRAPHIC_GMT

int gctl::regular_grid::get_xdim() const
{
	if (!initialized)
	{
		throw std::runtime_error("[gctl::regular_grid] The mesh is not initialized.");
	}

	return rg_xnum;
}

int gctl::regular_grid::get_ydim() const
{
	if (!initialized)
	{
		throw std::runtime_error("[gctl::regular_grid] The mesh is not initialized.");
	}

	return rg_ynum;
}

double gctl::regular_grid::get_xmin() const
{
	if (!initialized)
	{
		throw std::runtime_error("[gctl::regular_grid] The mesh is not initialized.");
	}

	return rg_xmin;
}

double gctl::regular_grid::get_ymin() const
{
	if (!initialized)
	{
		throw std::runtime_error("[gctl::regular_grid] The mesh is not initialized.");
	}

	return rg_ymin;
}

double gctl::regular_grid::get_dx() const
{
	if (!initialized)
	{
		throw std::runtime_error("[gctl::regular_grid] The mesh is not initialized.");
	}

	return rg_dx;
}

double gctl::regular_grid::get_dy() const
{
	if (!initialized)
	{
		throw std::runtime_error("[gctl::regular_grid] The mesh is not initialized.");
	}

	return rg_dy;
}

#ifdef GCTL_NETCDF

void gctl::regular_grid::load_netcdf_grid(std::string filename, mesh_data_type_e d_type, 
	std::string xname, std::string yname, std::string zname)
{
	gctl::array<double> in_x;
	gctl::array<double> in_y;
	gctl::array<std::string> in_name;
	gctl::matrix<double> in_arr;
	gctl::array<double> in_arr_single;

	// 如果未初始化 则初始化
	if (!initialized)
	{
		// read a data array
		read_netcdf_axis(filename, in_x, xname);
		read_netcdf_axis(filename, in_y, yname);
		// read single or block data
		if (zname == "null") read_netcdf_grid(filename, in_arr, in_name, xname, yname);
		else
		{
			read_netcdf_grid(filename, in_arr_single, xname, yname, zname);
			// copy data to matrix
			in_name.resize(1, zname);
			in_arr.resize(1, in_arr_single.size());
			for (size_t i = 0; i < in_arr.col_size(); i++)
			{
				in_arr[0][i] = in_arr_single[i];
			}
		}

		int xnum = in_x.size();
		int ynum = in_y.size();
		double xmin = in_x[0], ymin = in_y[0];
		double xmax = in_x[0], ymax = in_y[0];

		for (int i = 1; i < xnum; i++)
		{
			xmin = GCTL_MIN(xmin, in_x[i]);
			xmax = GCTL_MAX(xmax, in_x[i]);
		}

		for (int i = 1; i < ynum; i++)
		{
			ymin = GCTL_MIN(ymin, in_y[i]);
			ymax = GCTL_MAX(ymax, in_y[i]);
		}

		double dx = (xmax - xmin)/(xnum - 1);
		double dy = (ymax - ymin)/(ynum - 1);

		meshdata *data_ptr;
		if (d_type == NodeData)
		{
			init(filename, "Imported from a .nc file.", xnum, ynum, xmin, ymin, dx, dy);
			for (int i = 0; i < in_name.size(); i++)
			{
				if (!saved(in_name[i]))
				{
					data_ptr = meshdata_scalar::create(in_name[i], NodeData, xnum*ynum, true, 0.0);
					saved_data.push_back(data_ptr);
				}
				else throw std::runtime_error("[gctl::regular_grid] Duplicated data names.");
			}
		}
		else if (d_type == ElemData)
		{
			init(filename, "Imported from a .nc file.", xnum+1, ynum+1, xmin-0.5*dx, ymin-0.5*dy, dx, dy);
			for (int i = 0; i < in_name.size(); i++)
			{
				if (!saved(in_name[i]))
				{
					data_ptr = meshdata_scalar::create(in_name[i], ElemData, xnum*ynum, true, 0.0);
					saved_data.push_back(data_ptr);
				}
				else throw std::runtime_error("[gctl::regular_grid] Duplicated data names.");
			}
		}

		// get data pointer
		gctl::array<double> *val_ptr;
		for (size_t r = 0; r < in_arr.row_size(); r++)
		{
			val_ptr = (gctl::array<double>*) get_datval(in_name[r]);
			for (int i = 0; i < val_ptr->size(); i++)
			{
				val_ptr->at(i) = in_arr[r][i];
			}
		}

		initialized = true;
	}
	else
	{
		// read single or block data
		if (zname == "null") read_netcdf_grid(filename, in_arr, in_name, xname, yname);
		else
		{
			read_netcdf_grid(filename, in_arr_single, xname, yname, zname);
			// copy data to matrix
			in_name.resize(1, zname);
			in_arr.resize(1, in_arr_single.size());
			for (size_t i = 0; i < in_arr.col_size(); i++)
			{
				in_arr[0][i] = in_arr_single[i];
			}
		}

		meshdata *data_ptr;
		if (d_type == NodeData)
		{
			if (in_arr.col_size() != rg_xnum*rg_ynum)
			{
				throw std::runtime_error("[gctl::regular_grid] Incompatible input data size.");
			}

			for (int i = 0; i < in_name.size(); i++)
			{
				if (!saved(in_name[i]))
				{
					data_ptr = meshdata_scalar::create(in_name[i], d_type, rg_xnum*rg_ynum, true, 0.0);
					saved_data.push_back(data_ptr);
				}
				else throw std::runtime_error("[gctl::regular_grid] Duplicated data names.");
			}
		}
		else if (d_type == ElemData)
		{
			if (in_arr.col_size() != (rg_xnum-1)*(rg_ynum-1))
			{
				throw std::runtime_error("[gctl::regular_grid] Incompatible input data size.");
			}

			for (int i = 0; i < in_name.size(); i++)
			{
				if (!saved(in_name[i]))
				{
					data_ptr = meshdata_scalar::create(in_name[i], d_type, (rg_xnum-1)*(rg_ynum-1), true, 0.0);
					saved_data.push_back(data_ptr);
				}
				else throw std::runtime_error("[gctl::regular_grid] Duplicated data names.");
			}
		}

		// get data pointer
		gctl::array<double> *val_ptr;
		for (size_t r = 0; r < in_arr.row_size(); r++)
		{
			val_ptr = (gctl::array<double>*) get_datval(in_name[r]);
			for (int i = 0; i < val_ptr->size(); i++)
			{
				val_ptr->at(i) = in_arr[r][i];
			}
		}
	}
	return;
}

void gctl::regular_grid::save_netcdf_grid(std::string filename, mesh_data_type_e d_type, 
	std::string xname, std::string yname)
{
	if (!initialized)
	{
		throw std::runtime_error("[gctl::regular_grid] The mesh is not initialized.");
	}

	meshdata *curr_data = nullptr;

	bool if_append = false;
	for (iter = saved_data.begin(); iter != saved_data.end(); ++iter)
	{
		curr_data = *iter;
		if (curr_data->get_dattype() == d_type && d_type == NodeData && 
			curr_data->get_valtype() == Scalar && curr_data->get_output())
		{
			array<double>* data_ptr = (array<double>*)curr_data->get_datval_ptr();

			if (!if_append)
			{
				gctl::save_netcdf_grid(filename, *data_ptr, rg_xnum, rg_ynum, rg_xmin, 
					rg_dx, rg_ymin, rg_dy, xname, yname, curr_data->get_datname());
				if_append = true;
			}
			else
			{
				gctl::append_netcdf_grid(filename, *data_ptr, xname, yname, curr_data->get_datname());
			}
		}
		else if (curr_data->get_dattype() == d_type && d_type == ElemData && 
			curr_data->get_valtype() == Scalar && curr_data->get_output())
		{
			array<double>* data_ptr = (array<double>*)curr_data->get_datval_ptr();

			if (!if_append)
			{
				gctl::save_netcdf_grid(filename, *data_ptr, rg_xnum-1, rg_ynum-1, rg_xmin+0.5*rg_dx, 
					rg_dx, rg_ymin+0.5*rg_dy, rg_dy, xname, yname, curr_data->get_datname());
				if_append = true;
			}
			else
			{
				gctl::append_netcdf_grid(filename, *data_ptr, xname, yname, curr_data->get_datname());
			}
		}
	}

	return;
}

void gctl::regular_grid::save_netcdf_grid(std::string filename, std::string datname, 
	std::string xname, std::string yname)
{
	meshdata *curr_data = get_data(datname);

	if (!curr_data->get_output())
	{
		throw std::runtime_error("[gctl::regular_grid] Output of the data is disabled.");
	}

	if (curr_data->get_valtype() != Scalar)
	{
		throw std::runtime_error("[gctl::regular_grid] The output data type can't be saved to a netCDF file.");
	}

	if (curr_data->get_valtype() == Scalar && curr_data->get_dattype() == NodeData)
	{
		array<double>* data_ptr = (array<double>*)curr_data->get_datval_ptr();
		gctl::save_netcdf_grid(filename, *data_ptr, rg_xnum, rg_ynum, rg_xmin, 
			rg_dx, rg_ymin, rg_dy, xname, yname, curr_data->get_datname());
	}
	else if (curr_data->get_valtype() == Scalar && curr_data->get_dattype() == ElemData)
	{
		array<double>* data_ptr = (array<double>*)curr_data->get_datval_ptr();
		gctl::save_netcdf_grid(filename, *data_ptr, rg_xnum-1, rg_ynum-1, rg_xmin+0.5*rg_dx, 
			rg_dx, rg_ymin+0.5*rg_dy, rg_dy, xname, yname, curr_data->get_datname());
	}

	return;
}

#endif // GCTL_NETCDF

void gctl::regular_grid::load_surfer_grid(std::string filename, std::string datname, mesh_data_type_e d_type, surfer_file_type_e grid_type)
{
	gctl::array<double> in_arr;

	int xnum, ynum;
	double xmin, xmax, ymin, ymax, zmin, zmax, dx, dy, blank_val;

	// 如果微初始化 则初始化
	if (!initialized)
	{
		if (grid_type == Surfer6Text)
		{
			gctl::read_surfer6_grid(filename, in_arr, xnum, ynum, xmin, xmax, ymin, ymax, zmin, 
				zmax, gctl::Surfer6Text);

			dx = (xmax-xmin)/(xnum-1);
			dy = (ymax-ymin)/(ynum-1);
		}
		else if (grid_type == Surfer6Binary)
		{
			gctl::read_surfer6_grid(filename, in_arr, xnum, ynum, xmin, xmax, ymin, ymax, zmin, 
				zmax, gctl::Surfer6Binary);

			dx = (xmax-xmin)/(xnum-1);
			dy = (ymax-ymin)/(ynum-1);
		}
		else if (grid_type == Surfer7Grid)
		{
			gctl::read_surfer7_grid(filename, in_arr, xnum, ynum, xmin, ymin, dx, dy, zmin, zmax, blank_val);
		}

		meshdata *data_ptr;
		if (d_type == NodeData)
		{
			init(filename, "Imported from a surfer grid.", xnum, ynum, xmin, ymin, dx, dy);
			data_ptr = meshdata_scalar::create(datname, NodeData, xnum*ynum, true, 0.0);
		}
		else if (d_type == ElemData)
		{
			init(filename, "Imported from a surfer grid.", xnum+1, ynum+1, xmin-0.5*dx, ymin-0.5*dy, dx, dy);
			data_ptr = meshdata_scalar::create(datname, ElemData, xnum*ynum, true, 0.0);
		}

		saved_data.push_back(data_ptr);

		// get data pointer
		gctl::array<double> *val_ptr = (gctl::array<double>*) data_ptr->get_datval_ptr();
		for (int i = 0; i < val_ptr->size(); i++)
		{
			val_ptr->at(i) = in_arr[i];
		}

		initialized = true;
	}
	else
	{
		if (grid_type == Surfer6Text)
		{
			gctl::read_surfer6_grid(filename, in_arr, xnum, ynum, xmin, xmax, ymin, ymax, zmin, 
				zmax, gctl::Surfer6Text);
		}
		else if (grid_type == Surfer6Binary)
		{
			gctl::read_surfer6_grid(filename, in_arr, xnum, ynum, xmin, xmax, ymin, ymax, zmin, 
				zmax, gctl::Surfer6Binary);
		}
		else if (grid_type == Surfer7Grid)
		{
			gctl::read_surfer7_grid(filename, in_arr, xnum, ynum, xmin, ymin, dx, dy, zmin, zmax, blank_val);
		}

		meshdata *data_ptr;
		if (d_type == NodeData)
		{
			if (in_arr.size() != rg_xnum*rg_ynum)
			{
				throw std::runtime_error("[gctl::regular_grid] Incompatible input data size.");
			}

			data_ptr = meshdata_scalar::create(datname, d_type, rg_xnum*rg_ynum, true, 0.0);
		}
		else if (d_type == ElemData)
		{
			if (in_arr.size() != (rg_xnum-1)*(rg_ynum-1))
			{
				throw std::runtime_error("[gctl::regular_grid] Incompatible input data size.");
			}

			data_ptr = meshdata_scalar::create(datname, d_type, (rg_xnum-1)*(rg_ynum-1), true, 0.0);
		}

		saved_data.push_back(data_ptr);

		// get data pointer
		gctl::array<double> *val_ptr = (gctl::array<double>*) data_ptr->get_datval_ptr();
		for (int i = 0; i < val_ptr->size(); i++)
		{
			val_ptr->at(i) = in_arr[i];
		}
	}
	return;
}

void gctl::regular_grid::save_surfer_grid(std::string filename, std::string datname, surfer_file_type_e grid_type)
{
	meshdata *curr_data = curr_data = get_data(datname);

	if (!curr_data->get_output())
	{
		throw std::runtime_error("[gctl::regular_grid] Output of the data is disabled.");
	}

	if (curr_data->get_valtype() != Scalar)
	{
		throw std::runtime_error("[gctl::regular_grid] The output data type can't be saved to a Surfer file.");
	}

	gctl::array<double> *data_ptr = (gctl::array<double>*) curr_data->get_datval_ptr();
	if (curr_data->get_dattype() == NodeData)
	{
		if (grid_type == Surfer6Text)
		{
			gctl::save_surfer6_grid(filename, *data_ptr, rg_xnum, rg_ynum, 
				rg_xmin, rg_xmin+(rg_xnum-1)*rg_dx, rg_ymin, rg_ymin+(rg_ynum-1)*rg_dy, 
				NAN, NAN, gctl::Surfer6Text);
		}
		else if (grid_type == Surfer6Binary)
		{
			gctl::save_surfer6_grid(filename, *data_ptr, rg_xnum, rg_ynum, 
				rg_xmin, rg_xmin+(rg_xnum-1)*rg_dx, rg_ymin, rg_ymin+(rg_ynum-1)*rg_dy, 
				NAN, NAN, gctl::Surfer6Binary);
		}
		else if (grid_type == Surfer7Grid)
		{
			gctl::save_surfer6_grid(filename, *data_ptr, rg_xnum, rg_ynum, 
				rg_xmin, rg_ymin, rg_dx, rg_dy);
		}
	}
	else if (curr_data->get_dattype() == ElemData)
	{
		if (grid_type == Surfer6Text)
		{
			gctl::save_surfer6_grid(filename, *data_ptr, rg_xnum-1, rg_ynum-1, 
				rg_xmin+0.5*rg_dx, rg_xmin+0.5*rg_dx+(rg_xnum-2)*rg_dx, 
				rg_ymin+0.5*rg_dy, rg_ymin+0.5*rg_dy+(rg_ynum-2)*rg_dy, 
				NAN, NAN, gctl::Surfer6Text);
		}
		else if (grid_type == Surfer6Binary)
		{
			gctl::save_surfer6_grid(filename, *data_ptr, rg_xnum-1, rg_ynum-1, 
				rg_xmin+0.5*rg_dx, rg_xmin+0.5*rg_dx+(rg_xnum-2)*rg_dx, 
				rg_ymin+0.5*rg_dy, rg_ymin+0.5*rg_dy+(rg_ynum-2)*rg_dy, 
				NAN, NAN, gctl::Surfer6Binary);
		}
		else if (grid_type == Surfer7Grid)
		{
			gctl::save_surfer6_grid(filename, *data_ptr, rg_xnum-1, rg_ynum-1, 
				rg_xmin+0.5*rg_dx, rg_ymin+0.5*rg_dy, rg_dx, rg_dy);
		}
	}

	return;
}

void gctl::regular_grid::save_gmsh(std::string filename, index_packed_e packed)
{
	std::ofstream outfile;
	gctl::open_outfile(outfile, filename, ".msh");
	gctl::save2gmsh(outfile, elements, nodes, packed);
	outfile.close();
	return;
}

void gctl::regular_grid::save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed)
{
	base_mesh::save_gmsh(filename, d_type, out_mode, packed);
	return;
}

void gctl::regular_grid::save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed)
{
	base_mesh::save_gmsh(filename, datname, out_mode, packed);
	return;
}

void gctl::regular_grid::save_text(std::string filename, const array<std::string> &datname)
{
	mesh_data_type_e d_type;
	meshdata *dat_ptr = nullptr;

	dat_ptr = get_data(datname[0]);
	d_type = dat_ptr->get_dattype();
	for (size_t i = 1; i < datname.size(); i++)
	{
		dat_ptr = get_data(datname[i]);
		if (d_type != dat_ptr->get_dattype())
		{
			throw std::runtime_error("[gctl::regular_grid] multiple data location types occured.");
		}
	}

	array<array<double>*> datval_ptr(datname.size(), nullptr);
	for (size_t i = 0; i < datname.size(); i++)
	{
		datval_ptr[i] = (array<double>*) get_datval(datname[i]);
	}

	std::ofstream ofile;
	open_outfile(ofile, filename, ".txt");

	ofile << "# x y";
	for (size_t i = 0; i < datname.size(); i++)
	{
		ofile << " " << datname[i];
	}
	ofile << "\n";

	if (dat_ptr->get_dattype() == NodeData)
	{
		for (size_t i = 0; i < node_num; i++)
		{
			ofile << nodes[i].x << " " << nodes[i].y;
			for (size_t n = 0; n < datname.size(); n++)
			{
				ofile << " " << datval_ptr[n]->at(i);
			}
			ofile << "\n";
		}
	}
	else
	{
		for (size_t i = 0; i < ele_num; i++)
		{
			ofile << 0.5*(elements[i].dl->x + elements[i].ur->x) << " " << 0.5*(elements[i].dl->y + elements[i].ur->y);
			for (size_t n = 0; n < datname.size(); n++)
			{
				ofile << " " << datval_ptr[n]->at(i);
			}
			ofile << "\n";
		}
	}

	ofile.close();
	return;
}

void gctl::regular_grid::load_data_cloud(const array<point2dc> &in_posi, const array<double> &in_val, 
	double search_xlen, double search_ylen, double search_deg, std::string datname, mesh_data_type_e d_type)
{
	if (!initialized)
	{
		throw std::runtime_error("[gctl::regular_grid] The mesh is not initialized.");
	}

	gctl::array<gctl::point2dc> grid_points;

	if (in_posi.size() != in_val.size())
	{
		throw std::runtime_error("[gctl::regular_grid] Invalid sizes of the input data cloud.");
	}
	// create a new data object
	meshdata *data_ptr = add_data(datname, d_type, true, Scalar);

	array<double> *inte_data = nullptr;
	if (d_type == NodeData)
	{
		grid_points.resize(rg_xnum*rg_ynum);
		for (int j = 0; j < rg_ynum; j++)
		{
			for (int i = 0; i < rg_xnum; i++)
			{
				grid_points[i+j*rg_xnum].x = rg_xmin+i*rg_dx;
				grid_points[i+j*rg_xnum].y = rg_ymin+j*rg_dy;
			}
		}
	}
	else if (d_type == ElemData)
	{
		grid_points.resize((rg_xnum-1)*(rg_ynum-1));
		for (int j = 0; j < rg_ynum-1; j++)
		{
			for (int i = 0; i < rg_xnum-1; i++)
			{
				grid_points[i+j*(rg_xnum-1)].x = rg_xmin+i*rg_dx+0.5*rg_dx;
				grid_points[i+j*(rg_xnum-1)].y = rg_ymin+j*rg_dy+0.5*rg_dy;
			}
		}
	}

	inte_data = p2p_dist_2d(grid_points.get(), grid_points.size(), in_posi.get(), 
		in_val.get(), in_posi.size(), search_xlen, search_ylen, search_deg);

	array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
	for (int i = 0; i < val_ptr->size(); i++)
	{
		val_ptr->at(i) = inte_data->at(i);
	}

	delete inte_data;
	return;
}

void gctl::regular_grid::extract_points(std::string datname, const array<point2dc> &in_posi, 
	array<double> &out_val)
{
	if (in_posi.size() == 0)
	{
		throw std::runtime_error("[gctl::regular_grid] The input array is empty.");
	}

	meshdata *data_ptr = get_data(datname);
	mesh_data_value_e value_type = data_ptr->get_valtype();
	mesh_data_type_e data_type = data_ptr->get_dattype();
	if(value_type != Scalar)
	{
		throw std::runtime_error("[gctl::regular_grid] Wrong data's value type.");
	}

	array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
	out_val.resize(in_posi.size());

	int tmp_i, tmp_j, in_xnum, in_ynum;
	if (data_type == NodeData)
	{
		in_xnum = rg_xnum;
		in_ynum = rg_ynum;
		for (int i = 0; i < in_posi.size(); i++)
		{
			tmp_i = floor((in_posi[i].x - rg_xmin)/rg_dx);
			tmp_j = floor((in_posi[i].y - rg_ymin)/rg_dy);

			if (tmp_i > in_xnum-2 || tmp_i < 0 || tmp_j > in_ynum-2 || tmp_j < 0)
			{
				out_val[i] = NAN;
			}

			if (!std::isnan(val_ptr->at(tmp_i + tmp_j*in_xnum)) && 
				!std::isnan(val_ptr->at(tmp_i + 1 + tmp_j*in_xnum)) && 
				!std::isnan(val_ptr->at(tmp_i + 1 + (tmp_j+1)*in_xnum)) && 
				!std::isnan(val_ptr->at(tmp_i + (tmp_j+1)*in_xnum)))
			{
				out_val[i] = rect_interpolate(rg_xmin + tmp_i*rg_dx, rg_ymin + tmp_j*rg_dy, 
					rg_dx, rg_dy, in_posi[i].x, in_posi[i].y, val_ptr->at(tmp_i + tmp_j*in_xnum), 
					val_ptr->at(tmp_i + 1 + tmp_j*in_xnum), val_ptr->at(tmp_i + 1 + (tmp_j+1)*in_xnum), 
					val_ptr->at(tmp_i + (tmp_j+1)*in_xnum));
			}
			else out_val[i] = NAN;
		}
	}
	else
	{
		in_xnum = rg_xnum-1;
		in_ynum = rg_ynum-1;
		for (int i = 0; i < in_posi.size(); i++)
		{
			tmp_i = floor((in_posi[i].x - rg_xmin - 0.5*rg_dx)/rg_dx);
			tmp_j = floor((in_posi[i].y - rg_ymin - 0.5*rg_dy)/rg_dy);

			if (tmp_i > in_xnum-2 || tmp_i < 0 || tmp_j > in_ynum-2 || tmp_j < 0)
			{
				out_val[i] = NAN;
			}

			if (!std::isnan(val_ptr->at(tmp_i + tmp_j*in_xnum)) && 
				!std::isnan(val_ptr->at(tmp_i + 1 + tmp_j*in_xnum)) && 
				!std::isnan(val_ptr->at(tmp_i + 1 + (tmp_j+1)*in_xnum)) && 
				!std::isnan(val_ptr->at(tmp_i + (tmp_j+1)*in_xnum)))
			{
				out_val[i] = rect_interpolate(rg_xmin + 0.5*rg_dx + tmp_i*rg_dx, rg_ymin + 0.5*rg_dy + tmp_j*rg_dy, 
					rg_dx, rg_dy, in_posi[i].x, in_posi[i].y, val_ptr->at(tmp_i + tmp_j*in_xnum), 
					val_ptr->at(tmp_i + 1 + tmp_j*in_xnum), val_ptr->at(tmp_i + 1 + (tmp_j+1)*in_xnum), 
					val_ptr->at(tmp_i + (tmp_j+1)*in_xnum));
			}
			else out_val[i] = NAN;
		}
	}
	return;
}

void gctl::regular_grid::extract_profile(std::string datname, const point2dc &start_p, 
	const point2dc &end_p, int size_p, array<point2dc> &out_posi, array<double> &out_val)
{
	if (start_p == end_p)
	{
		throw std::runtime_error("[gctl::regular_grid] The start and end points can's be the same.");
	}

	if (size_p <= 0)
	{
		throw std::runtime_error("[gctl::regular_grid] Invalid profile size.");
	}

	point2dc dp = 1.0/(size_p - 1)*(end_p - start_p);
	out_posi.resize(size_p);
	for (size_t i = 0; i < size_p; i++)
	{
		out_posi[i] = start_p + (1.0*i)*dp;
	}

	extract_points(datname, out_posi, out_val);
	return;
}

void gctl::regular_grid::gradient(std::string datname, std::string gradname, gradient_type_e d_type, int order)
{
	meshdata *data_ptr = get_data(datname);
	mesh_data_value_e value_type = data_ptr->get_valtype();
	mesh_data_type_e data_type = data_ptr->get_dattype();

	if(value_type != Scalar)
	{
		throw std::runtime_error("[gctl::regular_grid] Wrong data's value type.");
	}

	if (d_type != Dx && d_type != Dy)
	{
		throw std::runtime_error("[gctl::regular_grid] The calculation type must be Dx or Dy.");
	}

	// 检查是否存在与梯度数据同名的数据 若有则检查是否需要重新建立数据
	meshdata *new_data_ptr = nullptr;
	mesh_data_value_e value_type2;
	mesh_data_type_e data_type2;
	if (saved(gradname))
	{
		new_data_ptr = get_data(gradname);
		value_type2 = new_data_ptr->get_valtype();
		data_type2  = new_data_ptr->get_dattype();
		if (value_type2 != value_type || data_type2 != data_type)
		{
			// 删除原有数据
			remove_data(gradname);
			// 新建数据
			new_data_ptr = add_data(gradname, data_type, true, value_type);
		}
	}
	else new_data_ptr = add_data(gradname, data_type, true, value_type);

	array<double> *inval_ptr = (array<double>*) data_ptr->get_datval_ptr();
	array<double> *outval_ptr = (array<double>*) new_data_ptr->get_datval_ptr();
	if (data_type == NodeData && d_type == Dx)
	{
		difference_2d(*inval_ptr, *outval_ptr, rg_ynum, rg_xnum, rg_dx, Dx, order);
	}
	else if (data_type == NodeData && d_type == Dy)
	{
		difference_2d(*inval_ptr, *outval_ptr, rg_ynum, rg_xnum, rg_dy, Dy, order);
	}
	else if (data_type == ElemData && d_type == Dx)
	{
		difference_2d(*inval_ptr, *outval_ptr, rg_ynum-1, rg_xnum-1, rg_dx, Dx, order);
	}
	else
	{
		difference_2d(*inval_ptr, *outval_ptr, rg_ynum-1, rg_xnum-1, rg_dy, Dy, order);
	}
	return;
}

void gctl::regular_grid::sum(std::string newname, std::string datname, std::string datname2)
{
	meshdata *data_ptr1 = get_data(datname);
	mesh_data_value_e value_type1 = data_ptr1->get_valtype();
	mesh_data_type_e data_type1 = data_ptr1->get_dattype();

	meshdata *data_ptr2 = get_data(datname2);
	mesh_data_value_e value_type2 = data_ptr2->get_valtype();
	mesh_data_type_e data_type2 = data_ptr2->get_dattype();

	if (value_type1 != value_type2 || data_type1 != data_type2)
	{
		throw std::runtime_error("[gctl::regular_grid] Can't preform the process.");
	}

	array<double> *val_ptr1 = (array<double>*) data_ptr1->get_datval_ptr();
	array<double> *val_ptr2 = (array<double>*) data_ptr2->get_datval_ptr();
	
	meshdata *new_data_ptr = nullptr;
	new_data_ptr = add_data(newname, data_type1, true, value_type1);

	array<double> *val_ptr = (array<double>*) new_data_ptr->get_datval_ptr();
	for (size_t i = 0; i < val_ptr->size(); i++)
	{
		val_ptr->at(i) = val_ptr1->at(i) + val_ptr2->at(i);
	}
	return;
}

void gctl::regular_grid::diff(std::string newname, std::string datname, std::string datname2)
{
	meshdata *data_ptr1 = get_data(datname);
	mesh_data_value_e value_type1 = data_ptr1->get_valtype();
	mesh_data_type_e data_type1 = data_ptr1->get_dattype();

	meshdata *data_ptr2 = get_data(datname2);
	mesh_data_value_e value_type2 = data_ptr2->get_valtype();
	mesh_data_type_e data_type2 = data_ptr2->get_dattype();

	if (value_type1 != value_type2 || data_type1 != data_type2)
	{
		throw std::runtime_error("[gctl::regular_grid] Can't preform the process.");
	}

	array<double> *val_ptr1 = (array<double>*) data_ptr1->get_datval_ptr();
	array<double> *val_ptr2 = (array<double>*) data_ptr2->get_datval_ptr();
	
	meshdata *new_data_ptr = nullptr;
	new_data_ptr = add_data(newname, data_type1, true, value_type1);

	array<double> *val_ptr = (array<double>*) new_data_ptr->get_datval_ptr();
	for (size_t i = 0; i < val_ptr->size(); i++)
	{
		val_ptr->at(i) = val_ptr1->at(i) - val_ptr2->at(i);
	}
	return;
}

void gctl::regular_grid::boolean(std::string newname, std::string datname, std::string maskname, bool reverse)
{
	meshdata *data_ptr1 = get_data(datname);
	mesh_data_value_e value_type1 = data_ptr1->get_valtype();
	mesh_data_type_e data_type1 = data_ptr1->get_dattype();

	meshdata *data_ptr2 = get_data(maskname);
	mesh_data_value_e value_type2 = data_ptr2->get_valtype();
	mesh_data_type_e data_type2 = data_ptr2->get_dattype();

	if (value_type1 != value_type2 || data_type1 != data_type2)
	{
		throw std::runtime_error("[gctl::regular_grid] Can't preform the process.");
	}

	array<double> *val_ptr1 = (array<double>*) data_ptr1->get_datval_ptr();
	array<double> *val_ptr2 = (array<double>*) data_ptr2->get_datval_ptr();
	
	meshdata *new_data_ptr = nullptr;
	new_data_ptr = add_data(newname, data_type1, true, value_type1);

	array<double> *val_ptr = (array<double>*) new_data_ptr->get_datval_ptr();
	if (!reverse)
	{
		for (size_t i = 0; i < val_ptr->size(); i++)
		{
			if (fabs(val_ptr2->at(i)) < 1e-8) val_ptr->at(i) = GCTL_BDL_MAX;
			else val_ptr->at(i) = val_ptr1->at(i);
		}
	}
	else
	{
		for (size_t i = 0; i < val_ptr->size(); i++)
		{
			if (fabs(val_ptr2->at(i)) < 1e-8) val_ptr->at(i) = val_ptr1->at(i);
			else val_ptr->at(i) = GCTL_BDL_MAX;
		}
	}
	return;
}

#ifdef GCTL_MESH_EXPRTK

void gctl::regular_grid::function(std::string expression_str, std::string newname, std::string x_str, std::string y_str, std::string f_str)
{
	exprtk::symbol_table<double> symbol_table;
	// We have three variables here.
	// The first one is 'f' which represents the resultent data
	// The last two are 'x' and 'y' by default 
	array<double> var(3);
	symbol_table.add_variable(f_str, var[0]);
	symbol_table.add_variable(x_str, var[1]);
	symbol_table.add_variable(y_str, var[2]);

	exprtk::expression<double> expression;
	expression.register_symbol_table(symbol_table);

	exprtk::parser<double> parser;
	if (!parser.compile(expression_str, expression)) throw std::runtime_error("[gctl::function] Fail to compile the math expression.");

	meshdata* data_ptr = nullptr;
	
	if (!saved(newname)) data_ptr = add_data(newname, NodeData, true, Scalar);
	else data_ptr = get_data(newname);

	array<double>* val_ptr = (array<double>*) data_ptr->get_datval_ptr();

	if (data_ptr->get_dattype() == NodeData)
	{
		for (size_t i = 0; i < val_ptr->size(); i++)
		{
			var[0] = val_ptr->at(i);
			var[1] = nodes[i].x;
			var[2] = nodes[i].y;
			val_ptr->at(i) = expression.value();
		}
	}
	else
	{
		for (size_t i = 0; i < val_ptr->size(); i++)
		{
			var[0] = val_ptr->at(i);
			var[1] = 0.5*(elements[i].dl->x + elements[i].ur->x);
			var[2] = 0.5*(elements[i].dl->y + elements[i].ur->y);
			val_ptr->at(i) = expression.value();
		}
	}
	return;
}

void gctl::regular_grid::calculator(std::string expression_str, const array<std::string> &var_list, const array<std::string> &data_list)
{
	if (var_list.size() != data_list.size() || var_list.size() < 2) throw std::runtime_error("[gctl::regular_grid] Invalid array sizes.");

	exprtk::symbol_table<double> symbol_table;
	array<double> var(var_list.size());

	for (size_t i = 0; i < var.size(); i++)
	{
		symbol_table.add_variable(var_list[i], var[i]);
	}

	exprtk::expression<double> expression;
	expression.register_symbol_table(symbol_table);

	exprtk::parser<double> parser;
	if (!parser.compile(expression_str, expression)) throw std::runtime_error("[gctl::regular_grid] Fail to compile the math expression.");

	array<meshdata*> data_ptrs(var.size(), nullptr);
	mesh_data_value_e val_type1, val_type2;
	mesh_data_type_e  data_type1, data_type2;
	
	// Remeber we put the output at the first place.
	for (size_t i = 1; i < var.size(); i++)
	{
		data_ptrs[i] = get_data(data_list[i]);
	}
	
	val_type1 = data_ptrs[1]->get_valtype();
	data_type1= data_ptrs[1]->get_dattype();
	for (size_t i = 2; i < var.size(); i++)
	{
		val_type2 = data_ptrs[i]->get_valtype();
		data_type2= data_ptrs[i]->get_dattype();
		if (val_type1 != val_type2 || data_type1 != data_type2)
		{
			throw std::runtime_error("[gctl::regular_grid] Can't preform the calculation on selected data.");
		}
	}

	
	if (!saved(data_list[0])) data_ptrs[0] = add_data(data_list[0], data_type1, true, val_type1);
	else
	{
		data_ptrs[0] = get_data(data_list[0]);
		val_type2 = data_ptrs[0]->get_valtype();
		data_type2= data_ptrs[0]->get_dattype();
		if (val_type1 != val_type2 || data_type1 != data_type2)
		{
			throw std::runtime_error("[gctl::regular_grid] Can't preform the calculation on selected data.");
		}
	}

	array<array<double>*> val_ptrs(var.size());
	for (size_t i = 0; i < var.size(); i++)
	{
		val_ptrs[i] = (array<double>*) data_ptrs[i]->get_datval_ptr();
	}

	for (size_t i = 0; i < val_ptrs[0]->size(); i++)
	{
		for (size_t v = 0; v < var.size(); v++)
		{
			var[v] = val_ptrs[v]->at(i);
		}

		val_ptrs[0]->at(i) = expression.value();
	}
	return;
}

#endif // GCTL_MESH_EXPRTK

#ifdef GCTL_MESH_WAVELIB

void gctl::regular_grid::wavelet(std::string datname, std::string wavename, int order, bool summary)
{
	meshdata *data_ptr = get_data(datname);
	mesh_data_value_e value_type = data_ptr->get_valtype();
	mesh_data_type_e data_type = data_ptr->get_dattype();

	if(value_type != Scalar)
	{
		throw std::runtime_error("[gctl::regular_grid] Wrong data's value type.");
	}

	if (order <= 0)
	{
		throw std::runtime_error("[gctl::regular_grid] The decomposition level must be greater than zero.");
	}

	array<double> *inval_ptr = (array<double>*) data_ptr->get_datval_ptr();

	wave_object obj = nullptr;
	wt2_object wt = nullptr;
	obj = wave_init(wavename.c_str());
	if (data_type == NodeData) wt = wt2_init(obj, "dwt", rg_ynum, rg_xnum, order);
	else wt = wt2_init(obj, "dwt", rg_ynum-1, rg_xnum-1, order);

	double *wavecoeffs = nullptr;
	wavecoeffs = dwt2(wt, inval_ptr->get());
	gctl::array<double> coeff_copy(wt->outlength);

	for (int i = 0; i < wt->outlength; i++)
	{
		coeff_copy[i] = wavecoeffs[i];
		wavecoeffs[i] = 0.0;
	}

	int ir = wt->dimensions[0], ic = wt->dimensions[1];
	for (int i = 0; i < ir*ic; i++)
	{
		wavecoeffs[i] = coeff_copy[i];
	}

	data_ptr = add_data(datname+"_A"+std::to_string(wt->J), data_type, true, value_type);
	array<double> *outval_ptr = (array<double>*) data_ptr->get_datval_ptr();
	idwt2(wt, wavecoeffs, outval_ptr->get());

	int coeff_idx;
	for (int j = 0; j < wt->J; j++)
	{
		for (int i = 0; i < wt->outlength; i++)
		{
			wavecoeffs[i] = 0.0;
		}

		ir = wt->dimensions[2*j];
		ic = wt->dimensions[2*j+1];
		coeff_idx = wt->coeffaccess[3*j+1];
		for (int i = 0; i < 3*ir*ic; i++)
		{
			wavecoeffs[i+coeff_idx] = coeff_copy[i+coeff_idx];
		}

		data_ptr = add_data(datname+"_D"+std::to_string(wt->J - j), data_type, true, value_type);
		outval_ptr = (array<double>*) data_ptr->get_datval_ptr();
		idwt2(wt, wavecoeffs, outval_ptr->get());
	}

	if (summary) wt2_summary(wt);

	if (obj != nullptr) wave_free(obj);
	if (wt != nullptr) wt2_free(wt);
	if (wavecoeffs != nullptr) free(wavecoeffs);
	return;
}

#endif // GCTL_MESH_WAVELIB