gctl_mesh/lib/mesh/regular_grid.cpp

/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * 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)
{
	check_initiated(true); // check to see if the mesh is already initiated
	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 << "Node num: " << node_num_ << std::endl;
	os << "Elem num: " << ele_num_ << std::endl;
	os << "X-range: " << rg_xmin << " -> " << rg_xmin + (rg_xnum - 1)*rg_dx << "\n";
	os << "Y-range: " << rg_ymin << " -> " << rg_ymin + (rg_ynum - 1)*rg_dy << "\n";
	os << "Interval: " << rg_dx << ", " << rg_dy << "\n";
	os << "Dimension: " << rg_xnum << ", " << rg_ynum << "\n";
	return;
}

void gctl::regular_grid::load_binary(std::string filename)
{
	check_initiated(true); // check to see if the mesh is already initiated

	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)
{
	check_initiated(); // check to see if the mesh is not initiated

	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)
{
	meshdata &data = get_data(datname);

	mathgl_dens plt; ///< mathgl绘图对象
	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(data.datval_, 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)
{
	meshdata &data = get_data(datname);

	gmt_JX_single pic;
	pic.set_command("psconvert", "-A -TG -E300");
	pic.plot(datname, data.datval_, 
		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
{
	check_initiated(); // check to see if the mesh is not initiated
	return rg_xnum;
}

int gctl::regular_grid::get_ydim() const
{
	check_initiated(); // check to see if the mesh is not initiated
	return rg_ynum;
}

double gctl::regular_grid::get_xmin() const
{
	check_initiated(); // check to see if the mesh is not initiated
	return rg_xmin;
}

double gctl::regular_grid::get_ymin() const
{
	check_initiated(); // check to see if the mesh is not initiated
	return rg_ymin;
}

double gctl::regular_grid::get_dx() const
{
	check_initiated(); // check to see if the mesh is not initiated
	return rg_dx;
}

double gctl::regular_grid::get_dy() const
{
	check_initiated(); // check to see if the mesh is not initiated
	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 new_data;
		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]))
				{
					new_data.clear();
					new_data.create(NodeData, Scalar, node_num_, in_name[i], true, GCTL_BDL_MAX);
					for (size_t j = 0; j < node_num_; j++)
					{
						new_data.datval_[j] = in_arr[i][j];
					}
					
					datalist_.push_back(new_data);
				}
				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]))
				{
					new_data.clear();
					new_data.create(ElemData, Scalar, ele_num_, in_name[i], true, GCTL_BDL_MAX);
					for (size_t j = 0; j < ele_num_; j++)
					{
						new_data.datval_[j] = in_arr[i][j];
					}
					
					datalist_.push_back(new_data);
				}
				else throw std::runtime_error("[gctl::regular_grid] Duplicated data names.");
			}
		}

		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 new_data;
		if (d_type == NodeData)
		{
			for (int i = 0; i < in_name.size(); i++)
			{
				if (!saved(in_name[i]))
				{
					new_data.clear();
					new_data.create(NodeData, Scalar, node_num_, in_name[i], true, GCTL_BDL_MAX);
					for (size_t j = 0; j < node_num_; j++)
					{
						new_data.datval_[i] = in_arr[i][j];
					}
					
					datalist_.push_back(new_data);
				}
				else throw std::runtime_error("[gctl::regular_grid] Duplicated data names.");
			}
		}
		else if (d_type == ElemData)
		{
			for (int i = 0; i < in_name.size(); i++)
			{
				if (!saved(in_name[i]))
				{
					new_data.clear();
					new_data.create(ElemData, Scalar, ele_num_, in_name[i], true, GCTL_BDL_MAX);
					for (size_t j = 0; j < ele_num_; j++)
					{
						new_data.datval_[i] = in_arr[i][j];
					}
					
					datalist_.push_back(new_data);
				}
				else throw std::runtime_error("[gctl::regular_grid] Duplicated data names.");
			}
		}
	}
	return;
}

void gctl::regular_grid::save_netcdf_grid(std::string filename, mesh_data_type_e d_type, 
	std::string xname, std::string yname)
{
	check_initiated(); // check to see if the mesh is not initiated

	bool if_append = false;
	for (size_t i = 0; i < datalist_.size(); i++)
	{
		if (datalist_[i].loctype_ == d_type && d_type == NodeData && 
			datalist_[i].valtype_ == Scalar && datalist_[i].output_ok_)
		{
			if (!if_append)
			{
				gctl::save_netcdf_grid(filename, datalist_[i].datval_, rg_xnum, rg_ynum, rg_xmin, 
					rg_dx, rg_ymin, rg_dy, xname, yname, datalist_[i].name_);
				if_append = true;
			}
			else
			{
				gctl::append_netcdf_grid(filename, datalist_[i].datval_, xname, yname, datalist_[i].name_);
			}
		}
		else if (datalist_[i].loctype_ == d_type && d_type == ElemData && 
			datalist_[i].valtype_ == Scalar && datalist_[i].output_ok_)
		{
			if (!if_append)
			{
				gctl::save_netcdf_grid(filename, datalist_[i].datval_, rg_xnum-1, rg_ynum-1, rg_xmin+0.5*rg_dx, 
					rg_dx, rg_ymin+0.5*rg_dy, rg_dy, xname, yname, datalist_[i].name_);
				if_append = true;
			}
			else
			{
				gctl::append_netcdf_grid(filename, datalist_[i].datval_, xname, yname, datalist_[i].name_);
			}
		}
	}
	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.output_ok_)
	{
		throw std::runtime_error("[gctl::regular_grid] Output of the data is disabled.");
	}

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

	if (curr_data.valtype_ == Scalar && curr_data.loctype_ == NodeData)
	{
		gctl::save_netcdf_grid(filename, curr_data.datval_, rg_xnum, rg_ynum, rg_xmin, 
			rg_dx, rg_ymin, rg_dy, xname, yname, curr_data.name_);
	}
	else if (curr_data.valtype_ == Scalar && curr_data.loctype_ == ElemData)
	{
		gctl::save_netcdf_grid(filename, curr_data.datval_, 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.name_);
	}
	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)
{
	int xnum, ynum;
	double xmin, xmax, ymin, ymax, zmin, zmax, dx, dy, blank_val;
	gctl::array<double> in_arr;

	// 如果微初始化 则初始化
	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 curr_data;
		if (d_type == NodeData)
		{
			init(filename, "Imported from a surfer grid.", xnum, ynum, xmin, ymin, dx, dy);
			curr_data.create(NodeData, Scalar, node_num_, datname, true, GCTL_BDL_MAX);
			for (size_t i = 0; i < node_num_; i++)
			{
				curr_data.datval_[i] = in_arr[i];
			}
		}
		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);
			curr_data.create(ElemData, Scalar, ele_num_, datname, true, GCTL_BDL_MAX);
			for (size_t i = 0; i < ele_num_; i++)
			{
				curr_data.datval_[i] = in_arr[i];
			}
		}

		datalist_.push_back(curr_data);
		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 curr_data;
		if (d_type == NodeData)
		{
			if (in_arr.size() != rg_xnum*rg_ynum)
			{
				throw std::runtime_error("[gctl::regular_grid] Incompatible input data size.");
			}

			curr_data.create(NodeData, Scalar, node_num_, datname, true, GCTL_BDL_MAX);
			for (size_t i = 0; i < node_num_; i++)
			{
				curr_data.datval_[i] = in_arr[i];
			}
		}
		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.");
			}

			curr_data.create(ElemData, Scalar, ele_num_, datname, true, GCTL_BDL_MAX);
			for (size_t i = 0; i < ele_num_; i++)
			{
				curr_data.datval_[i] = in_arr[i];
			}
		}

		datalist_.push_back(curr_data);
	}
	return;
}

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

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

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

	if (curr_data.loctype_ == NodeData)
	{
		if (grid_type == Surfer6Text)
		{
			gctl::save_surfer6_grid(filename, curr_data.datval_, 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, curr_data.datval_, 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, curr_data.datval_, rg_xnum, rg_ynum, 
				rg_xmin, rg_ymin, rg_dx, rg_dy);
		}
	}
	else if (curr_data.loctype_ == ElemData)
	{
		if (grid_type == Surfer6Text)
		{
			gctl::save_surfer6_grid(filename, curr_data.datval_, 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, curr_data.datval_, 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, curr_data.datval_, 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, output_type_e out_mode, index_packed_e packed)
{
	base_mesh::save_gmsh_withdata(filename, 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_withdata(filename, datname, out_mode, packed);
	return;
}

void gctl::regular_grid::save_text(std::string filename, const array<std::string> &datname)
{
	meshdata &curr_data = get_data(datname[0]);
	for (size_t i = 1; i < datalist_.size(); i++)
	{
		if (curr_data.loctype_ != datalist_[i].loctype_)
		{
			throw std::runtime_error("[gctl::regular_grid] multiple data location types found.");
		}
	}

	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 (curr_data.loctype_ == NodeData)
	{
		for (size_t i = 0; i < node_num_; i++)
		{
			ofile << nodes[i].x << " " << nodes[i].y;
			for (size_t n = 0; n < datalist_.size(); n++)
			{
				ofile << " " << datalist_[n].datval_[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 < datalist_.size(); n++)
			{
				ofile << " " << datalist_[n].datval_[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)
{
	check_initiated(); // 检查网格是否已经初始化

	if (in_posi.size() != in_val.size())
	{
		throw std::runtime_error("[gctl::regular_grid] Invalid sizes of the input data cloud.");
	}

	array<double> *inte_data = nullptr;
	gctl::array<gctl::point2dc> grid_points;
	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);
	
	// create a new data object
	meshdata &curr_data = add_data(d_type, Scalar, datname, 0.0, true, GCTL_BDL_MAX);
	for (int i = 0; i < curr_data.datval_.size(); i++)
	{
		curr_data.datval_[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 &curr_data = get_data(datname);
	mesh_data_value_e value_type = curr_data.valtype_;
	mesh_data_type_e data_type = curr_data.loctype_;
	if(value_type != Scalar)
	{
		throw std::runtime_error("[gctl::regular_grid] Wrong data's value type.");
	}

	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(curr_data.datval_[tmp_i + tmp_j*in_xnum]) && 
				!std::isnan(curr_data.datval_[tmp_i + 1 + tmp_j*in_xnum]) && 
				!std::isnan(curr_data.datval_[tmp_i + 1 + (tmp_j+1)*in_xnum]) && 
				!std::isnan(curr_data.datval_[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, curr_data.datval_[tmp_i + tmp_j*in_xnum], 
					curr_data.datval_[tmp_i + 1 + tmp_j*in_xnum], curr_data.datval_[tmp_i + 1 + (tmp_j+1)*in_xnum], 
					curr_data.datval_[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(curr_data.datval_[tmp_i + tmp_j*in_xnum]) && 
				!std::isnan(curr_data.datval_[tmp_i + 1 + tmp_j*in_xnum]) && 
				!std::isnan(curr_data.datval_[tmp_i + 1 + (tmp_j+1)*in_xnum]) && 
				!std::isnan(curr_data.datval_[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, curr_data.datval_[tmp_i + tmp_j*in_xnum], 
					curr_data.datval_[tmp_i + 1 + tmp_j*in_xnum], curr_data.datval_[tmp_i + 1 + (tmp_j+1)*in_xnum], 
					curr_data.datval_[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 &curr_data = get_data(datname);
	mesh_data_value_e value_type = curr_data.valtype_;
	mesh_data_type_e data_type = curr_data.loctype_;

	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.");
	}

	// 检查是否存在与梯度数据同名的数据 若有则检查是否需要重新建立数据
	if (saved(gradname))
	{
		meshdata &check_data = get_data(gradname);
		mesh_data_value_e value_type2 = check_data.valtype_;
		mesh_data_type_e data_type2  = check_data.loctype_;
		// 删除原有数据
		if (value_type2 != value_type || data_type2 != data_type)
		{
			remove_data(gradname);
		}
		add_data(data_type, Scalar, gradname, 0.0, true, GCTL_BDL_MAX);
	}
	else add_data(data_type, Scalar, gradname, 0.0, true, GCTL_BDL_MAX);

	meshdata &new_data = get_data(gradname);
	if (data_type == NodeData && d_type == Dx)
	{
		difference_2d(curr_data.datval_, new_data.datval_, rg_ynum, rg_xnum, rg_dx, Dx, order);
	}
	else if (data_type == NodeData && d_type == Dy)
	{
		difference_2d(curr_data.datval_, new_data.datval_, rg_ynum, rg_xnum, rg_dy, Dy, order);
	}
	else if (data_type == ElemData && d_type == Dx)
	{
		difference_2d(curr_data.datval_, new_data.datval_, rg_ynum-1, rg_xnum-1, rg_dx, Dx, order);
	}
	else
	{
		difference_2d(curr_data.datval_, new_data.datval_, 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_1 = get_data(datname);
	mesh_data_value_e value_type1 = data_1.valtype_;
	mesh_data_type_e data_type1 = data_1.loctype_;

	meshdata &data_2 = get_data(datname2);
	mesh_data_value_e value_type2 = data_2.valtype_;
	mesh_data_type_e data_type2 = data_2.loctype_;

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

	meshdata &new_data = add_data(data_type1, value_type1, newname, 0.0, true, GCTL_BDL_MAX);

	for (size_t i = 0; i < new_data.datval_.size(); i++)
	{
		new_data.datval_[i] = data_1.datval_[i] + data_2.datval_[i];
	}
	return;
}

void gctl::regular_grid::diff(std::string newname, std::string datname, std::string datname2)
{
	meshdata &data_1 = get_data(datname);
	mesh_data_value_e value_type1 = data_1.valtype_;
	mesh_data_type_e data_type1 = data_1.loctype_;

	meshdata &data_2 = get_data(datname2);
	mesh_data_value_e value_type2 = data_2.valtype_;
	mesh_data_type_e data_type2 = data_2.loctype_;

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

	meshdata &new_data = add_data(data_type1, value_type1, newname, 0.0, true, GCTL_BDL_MAX);

	for (size_t i = 0; i < new_data.datval_.size(); i++)
	{
		new_data.datval_[i] = data_1.datval_[i] - data_2.datval_[i];
	}
	return;
}

void gctl::regular_grid::boolean(std::string newname, std::string datname, std::string maskname, bool reverse)
{
	meshdata &data_1 = get_data(datname);
	mesh_data_value_e value_type1 = data_1.valtype_;
	mesh_data_type_e data_type1 = data_1.loctype_;

	meshdata &data_2 = get_data(maskname);
	mesh_data_value_e value_type2 = data_2.valtype_;
	mesh_data_type_e data_type2 = data_2.loctype_;

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

	meshdata &new_data = add_data(data_type1, value_type1, newname, 0.0, true, GCTL_BDL_MAX);

	if (!reverse)
	{
		for (size_t i = 0; i < new_data.datval_.size(); i++)
		{
			if (fabs(data_2.datval_[i]) < 1e-8 || isnan(data_2.datval_[i])) new_data.datval_[i] = GCTL_BDL_MAX;
			else new_data.datval_[i] = data_1.datval_[i];
		}
	}
	else
	{
		for (size_t i = 0; i < new_data.datval_.size(); i++)
		{
			if (fabs(data_2.datval_[i]) < 1e-8 || isnan(data_2.datval_[i])) new_data.datval_[i] = data_1.datval_[i];
			else new_data.datval_[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.");

	if (!saved(newname)) add_data(NodeData, Scalar, newname, 0.0, true, GCTL_BDL_MAX);
	meshdata &curr_data = get_data(newname);

	if (curr_data.loctype_ == NodeData)
	{
		for (size_t i = 0; i < curr_data.datval_.size(); i++)
		{
			var[0] = curr_data.datval_[i];
			var[1] = nodes[i].x;
			var[2] = nodes[i].y;
			curr_data.datval_[i] = expression.value();
		}
	}
	else
	{
		for (size_t i = 0; i < curr_data.datval_.size(); i++)
		{
			var[0] = curr_data.datval_[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);
			curr_data.datval_[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*> datas_ptr(var.size());
	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++)
	{
		datas_ptr[i] = get_data_ptr(data_list[i]);
	}
	
	val_type1 = datas_ptr[1]->valtype_;
	data_type1= datas_ptr[1]->loctype_;
	for (size_t i = 2; i < var.size(); i++)
	{
		val_type2 = datas_ptr[i]->valtype_;
		data_type2= datas_ptr[i]->loctype_;
		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])) add_data(data_type1, val_type1, data_list[0], 0.0, true, GCTL_BDL_MAX);

	datas_ptr[0] = get_data_ptr(data_list[0]);
	val_type2 = datas_ptr[0]->valtype_;
	data_type2= datas_ptr[0]->loctype_;
	if (val_type1 != val_type2 || data_type1 != data_type2)
	{
		throw std::runtime_error("[gctl::regular_grid] Can't preform the calculation on selected data.");
	}

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

		datas_ptr[0]->datval_[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 &curr_data = get_data(datname);
	mesh_data_value_e value_type = curr_data.valtype_;
	mesh_data_type_e data_type = curr_data.loctype_;

	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.");
	}

	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, curr_data.datval_.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];
	}

	meshdata &new_data = add_data(data_type, value_type, datname+"_A"+std::to_string(wt->J), 0.0, true, GCTL_BDL_MAX);
	idwt2(wt, wavecoeffs, new_data.datval_.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];
		}

		new_data = add_data(data_type, value_type, datname+"_D"+std::to_string(wt->J - j), 0.0, true, GCTL_BDL_MAX);
		idwt2(wt, wavecoeffs, new_data.datval_.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