1290 lines
36 KiB
C++
1290 lines
36 KiB
C++
|
/********************************************************
|
|||
|
|
* ██████╗ ██████╗████████╗██╗
|
|||
|
|
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
|||
|
|
* ██║ ███╗██║ ██║ ██║
|
|||
|
|
* ██║ ██║██║ ██║ ██║
|
|||
|
|
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
|||
|
|
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
|||
|
|
* 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)
|
|||
|
|
{
|
|||
|
|
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, "x", "y", curr_data->get_datname());
|
|||
|
|
if_append = true;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
gctl::append_netcdf_grid(filename, *data_ptr, "x", "y", 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, "x", "y", curr_data->get_datname());
|
|||
|
|
if_append = true;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
gctl::append_netcdf_grid(filename, *data_ptr, "x", "y", curr_data->get_datname());
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
void gctl::regular_grid::save_netcdf_grid(std::string filename, std::string datname)
|
|||
|
|
{
|
|||
|
|
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, "x", "y", 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, "x", "y", 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.");
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
linespace(start_p, end_p, size_p, out_posi);
|
|||
|
|
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
|