gctl_toolkits/tetgm/tetgm_backup

/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2022  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.
 ******************************************************/

// add gctl head files
#include "gctl/io.h"
#include "gctl/potential.h"
#include "gctl/utility.h"

int main(int argc, char *argv[]) try
{
	const char* pro_name = "tetgm";
	const char* pro_info = "1.0 - Forward modeling of gravitational data using unstructured mesh. \
		This program is a toolkit of the GCTL package. The GCTL comes with ABSOLUTE NO WARRANTY. \
		Please see instructions or contact the author for more information.";
	const char* model_file_str = "Name of the 3D density model file. The program takes either Gmsh (.msh) or Tetgen model files. \
		For Tetgen files, one must have both .node and .ele files. NOTE: name extensions are not needed for Tetgen files.";
	const char* observation_str = "Name of the observation file. Or initializing parameters of the observation points. \
		The input file should have at least three data columns that represent x, y and z coordinates of the observation \
		points. Use -d option to select data columns (default is 0,1,2).";
	const char* physics_str = "Name of a file that contains the 3D model's physical properties. The file should at least has one column which \
		represents model densities. Use -d option to select data columns (default is 0). The option could also take name of the .msh \
		file's model data. Otherwise, entry a float value to indicate evenly distributed physical property.";
	const char* output_str = "Prefix of the output file's name. _Vz, _Vzx, _Vzy and _Vzz suffix will be added automatically \
		according to the gravitational component that is calculated.";

	char model_file[1024] = "NULL";
	char obser_para[1024] = "NULL";
	char physic_para[1024]= "NULL";
	char output_file[1024]= "NULL";
	gctl::index_packed_e compact_mode = gctl::NotPacked;
	bool components_sign[4] = {false, false, false, false};
	gctl::gravitational_field_type_e components_mark[4] = {gctl::Vz, gctl::Tzx, gctl::Tzy, gctl::Tzz};

	// option format
	// char* name, int has_arg, int *flag, int val, char* info, char* format, bool manda
	static struct gctl::option_info long_opts_help[] = 
	{
		{model_file_str, "<filename>", true},
		{observation_str, "<filename>[+d<x-col>,<y-col>,<z-col>]|<xmin>/<dx>/<xmax>/<ymin>/<dy>/<ymax>/<elevation>", true},
		{physics_str, "<filename>+d<col>|<dataname>|<value>", true},
		{output_str, "<filename>", true},
		{"Calculate gravity data.", 0, false},
		{"Calculate gravity gradient data along x-axis.", 0, false},
		{"Calculate gravity gradient data along y-axis.", 0, false},
		{"Calculate gravity gradient data along z-axis.", 0, false},
		{"Set this option if the index of the mesh's vertice and elements are starting from zero.", 0, false},
		{"Show help information.", 0, false},
		{0, 0, 0}
	};

	static struct option long_opts[] = 
	{
		{"model-file", required_argument, NULL, 'm'},
		{"observation", required_argument, NULL, 'b'},
		{"physics", required_argument, NULL, 'p'},
		{"output-file", required_argument, NULL, 'o'},
		{"gravity", no_argument, NULL, 'g'},
		{"gradient-x", no_argument, NULL, 'x'},
		{"gradient-y", no_argument, NULL, 'y'},
		{"gradient-z", no_argument, NULL, 'z'},
		{"compact", no_argument, NULL, 'c'},
		{"help", no_argument, NULL, 'h'},
		{0, 0, 0, 0}
	};

	int curr;
	while (1)
	{
		int opt_index = 0;

		if (argc == 1)
		{
			gctl::display_logo(std::clog);
			gctl::getopt_long_help(long_opts, long_opts_help, pro_name, pro_info);
			return 0;
		}

		curr = getopt_long(argc, argv, "gxyzchm:b:p:o:", long_opts, &opt_index);

		if (curr == -1) break;

		switch (curr)
		{
			case 'h': //显示帮助信息
				gctl::display_logo(std::clog);
				gctl::getopt_long_help(long_opts, long_opts_help, pro_name, pro_info);
				return 0;
			case 'm':
				sscanf(optarg, "%s", model_file);
				break;
			case 'b':
				sscanf(optarg, "%s", obser_para);
				break;
			case 'p':
				sscanf(optarg, "%s", physic_para);
				break;
			case 'o':
				sscanf(optarg, "%s", output_file);
				break;
			case 'g':
				components_sign[0] = true;
				break;
			case 'x':
				components_sign[1] = true;
				break;
			case 'y':
				components_sign[2] = true;
				break;
			case 'z':
				components_sign[3] = true;
				break;
			case 'c':
				compact_mode = gctl::Packed;
				break;
			case '?':
				gctl::getopt_long_help(long_opts, long_opts_help, pro_name, " - Unknown options. Please see the help information below.");
				break;
			default:
				abort();
		}
	}

	// check for mandatory option(s)
	if (!strcmp(model_file, "NULL") || !strcmp(obser_para, "NULL") || 
		!strcmp(physic_para, "NULL")|| !strcmp(output_file, "NULL"))
	{
		GCTL_ShowWhatError("Missing arguments for mandatory option(s).", GCTL_ERROR_ERROR, 
			0, "See tips below. Or use -h option to see the full instruction.", 0);
		if (!strcmp(model_file, "NULL"))
			gctl::getopt_long_option_info('m', long_opts, long_opts_help);
		if (!strcmp(obser_para, "NULL"))
			gctl::getopt_long_option_info('b', long_opts, long_opts_help);
		if (!strcmp(physic_para, "NULL"))
			gctl::getopt_long_option_info('p', long_opts, long_opts_help);
		if (!strcmp(output_file, "NULL"))
			gctl::getopt_long_option_info('o', long_opts, long_opts_help);
		return 0;
	}

	// check for forward modeling data type
	if (components_sign[0] == false && components_sign[1] == false && 
		components_sign[2] == false && components_sign[3] == false)
	{
		GCTL_ShowWhatError("No forward modeling data type is selected.", GCTL_ERROR_ERROR, 
			0, "Use -h option to see the full instruction.", 0);
		return 0;
	}

	// declare variables here
	bool gmsh_file = false;
	// we firstly try to read points
	gctl::array<gctl::point3dc> obs_points;
	// get tetrahedron mesh
	gctl::array<gctl::vertex3dc> mesh_node;
	gctl::array<gctl::grav_tetrahedron> mesh_tet;
	gctl::array<gctl::gravtet_para> mesh_tet_gp;
	// define physic array
	gctl::array<double> mesh_rho;
	// define result array
	gctl::array<gctl::point3dc> grav_grad;
	gctl::array<gctl::tensor> grav_tens;

	gctl::text_descriptor desc;

	// start the sequence here
	std::ifstream infile;
	// initialize observation points here
	double xmin, dx, xmax, ymin, dy, ymax, ele;
	if (7 == sscanf(obser_para, "%lf/%lf/%lf/%lf/%lf/%lf/%lf", 
		&xmin, &dx, &xmax, &ymin, &dy, &ymax, &ele))
	{
		gctl::get_grid_point3c(obs_points, xmin, xmax, ymin, ymax, dx, dy, ele);
	}
	else
	{
		char obs_file[1024], obs_order[1024];
		if (2 == sscanf(obser_para,"%[^+]+d%s", obs_file, obs_order))
		{
			gctl::get_xyz_points(obs_file, obs_points, desc, obs_order);
		}
		else
		{
			gctl::get_xyz_points(obser_para, obs_points, desc);
		}
	}

	std::string tmp_str = model_file;
	// proceed as Gmsh files
	if (tmp_str.substr(tmp_str.length()-4, tmp_str.length()) == ".msh")
	{
		gmsh_file = true;

		gctl::open_infile(infile, model_file);
		gctl::read_gmsh_node(infile, mesh_node, compact_mode);
		gctl::read_gmsh_element(infile, mesh_tet, mesh_node, compact_mode);
		infile.close();
	}
	// proceed as tetgen files
	else
	{
		gctl::read_Tetgen_node(model_file, mesh_node, compact_mode);
		gctl::read_Tetgen_element(model_file, mesh_tet, mesh_node, compact_mode);
	}

	// initiate gravtet_para
	gctl::callink_gravity_para(mesh_tet, mesh_tet_gp);

	// Firstly try to explain physic_para as float number
	double physic_value;
	if (1 == sscanf(physic_para, "%lf", &physic_value))
	{
		mesh_rho.resize(mesh_tet.size(), physic_value);
	}
	else
	{
		// try to use physic_para as file name
		std::vector<std::vector<double> > txt_content;
		try
		{
			char physic_filename[1024];
			//默认的读入的数据列为第一列
			int physic_col = 0;
			if (2 != sscanf(physic_para,"%[^+]+d%d", physic_filename, &physic_col))
				strcpy(physic_filename, physic_para);

			gctl::read_text2vector2d(physic_filename, txt_content, desc);

			if (txt_content.size() != mesh_tet.size())
			{
				std::string msg = "Element size doesn't match.";
				throw msg;
			}

			mesh_rho.resize(txt_content.size());
			for (int i = 0; i < mesh_rho.size(); i++)
				mesh_rho.at(i) = txt_content.at(i).at(physic_col);

			gctl::destroy_vector(txt_content);

			goto forward_calculation;
		}
		catch(std::string err_str)
		{
			if (!gmsh_file) throw err_str;
			else
			{
				char msg[1024] = "Using ";
				strcat(msg, physic_para);
				strcat(msg, " as a data name in ");
				strcat(msg, model_file);
				GCTL_ShowWhatError(err_str, GCTL_MESSAGE_ERROR, msg, 0, 0);
			}
		}
		// try to use physic_para as data_name
		std::ifstream infile;
		gctl::open_infile(infile, model_file);
		gctl::read_gmsh_data(infile, mesh_rho, physic_para);
		infile.close();
		if (mesh_rho.empty())
		{
			std::string msg = "Density model is not found!";
			throw msg;
		}
	}

forward_calculation:
	std::vector<std::vector<double> > save_content;
	save_content.resize(obs_points.size());
	for (int i = 0; i < obs_points.size(); i++)
		save_content[i].resize(4);

	for (int i = 0; i < save_content.size(); i++)
	{
		save_content[i][0] = obs_points.at(i).x;
		save_content[i][1] = obs_points.at(i).y;
		save_content[i][2] = obs_points.at(i).z;
	}

	std::string out_name = output_file;
	if (components_sign[0])
	{
		gctl::gobser(grav_grad, mesh_tet, obs_points, mesh_rho);
		for (int i = 0; i < obs_points.size(); ++i)
		{
			save_content[i][3] = grav_grad[i].z;
		}

		gctl::save_vector2d2text(out_name+"_Vz.txt", save_content);
	}

	if (components_sign[1] || components_sign[2] || components_sign[3])
	{
		gctl::gobser(grav_tens, mesh_tet, obs_points, mesh_rho);

		if (components_sign[1])
		{
			for (int i = 0; i < obs_points.size(); ++i)
			{
				save_content[i][3] = grav_tens[i].at(2, 0);
			}

			gctl::save_vector2d2text(out_name+"_Vzx.txt", save_content);
		}

		if (components_sign[2])
		{
			for (int i = 0; i < obs_points.size(); ++i)
			{
				save_content[i][3] = grav_tens[i].at(2, 1);
			}

			gctl::save_vector2d2text(out_name+"_Vzy.txt", save_content);
		}

		if (components_sign[3])
		{
			for (int i = 0; i < obs_points.size(); ++i)
			{
				save_content[i][3] = grav_tens[i].at(2, 2);
			}

			gctl::save_vector2d2text(out_name+"_Vzz.txt", save_content);
		}
	}
	return 0;
}
catch(std::exception &e)
{
	GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
}