gctl_toolkits/tetgm/tetgm.cpp

/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * 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/core.h"
#include "gctl/io.h"
#include "gctl/potential.h"
#include "gctl/utility.h"

using namespace gctl;

const char* pro_info_str = "1.1 - Forward modeling of the gravitational and magnetic data using \
tetrahedron meshs. 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 model file. Add .msh extension to specify a Gmsh model file. \
For Tetgen files, one must have both .node and .ele files. And name extensions are not needed for the 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.";

const char* mag_str = "Magnetization parameters. Needed for forward modeling magnetic data. The geo-magnetic field parameters are only used for total intense data.";

int main(int argc, char *argv[]) try
{
	flags_parser fp;
	fp.set_proname("tetgm");
	fp.set_proinfo(pro_info_str);
	fp.add_opt('m', "model-file", required_argument, NULL, model_file_str, "<file>", true);
	fp.add_opt('b', "observation", required_argument, NULL, observation_str, "<file>[+d<x-col>,<y-col>,<z-col>]|<xmin>/<dx>/<xmax>/<ymin>/<dy>/<ymax>/<elevation>", true);
	fp.add_opt('p', "physics", required_argument, NULL, physics_str, "<file>+d<col>|<Gmsh-data>|<value>", true);
	fp.add_opt('o', "output-file", required_argument, NULL, output_str, "<file>", true);
	fp.add_opt('a', "magnetization", required_argument, NULL, mag_str, "<mag-inc>/<mag-dec>/<geo-inc>/<geo-dec>", false);
	fp.add_opt('g', "gravity", no_argument, NULL, "Calculate gravity data.", 0, false);
	fp.add_opt('x', "gradient-x", no_argument, NULL, "Calculate gravity gradient data along x-axis.", 0, false);
	fp.add_opt('y', "gradient-y", no_argument, NULL, "Calculate gravity gradient data along y-axis.", 0, false);
	fp.add_opt('z', "gradient-z", no_argument, NULL, "Calculate gravity gradient data along z-axis.", 0, false);
	fp.add_opt('t', "magnetic", no_argument, NULL, "Calculate total intense magnetic data.", 0, false);
	fp.add_opt('u', "magnetic-x", no_argument, NULL, "Calculate the x-componment of the magnetic data.", 0, false);
	fp.add_opt('v', "magnetic-y", no_argument, NULL, "Calculate the y-componment of the magnetic data.", 0, false);
	fp.add_opt('r', "magnetic-z", no_argument, NULL, "Calculate the z-componment of the magnetic data.", 0, false);
	fp.add_opt('c', "compact", no_argument, NULL, "Set this option if the index of the mesh's vertice and elements are starting from zero.", 0, false);
	fp.add_opt('h', "help", no_argument, NULL, "Show help information.", 0, false);
	fp.configure(argc, argv);

	if (argc == 1 || fp.set_opt('h'))
	{
		fp.show_help_page();
		return 0;
	}

	std::string model_file, obser_para, physic_para, mag_para, output_file;
	gctl::index_packed_e compact_mode = gctl::NotPacked;
	bool components_sign[4] = {false, false, false, false};
	bool components_sign_mag[4] = {false, false, false, false};
	gctl::gravitational_field_type_e components_mark[4] = {gctl::Vz, gctl::Tzx, gctl::Tzy, gctl::Tzz};
	gctl::magnetic_field_type_e components_mark_mag[4] = {gctl::DeltaT, gctl::Bx, gctl::By, gctl::Bz};

	fp.get_argv(
		{'m', 'b', 'p', 'a', 'o'}, 
		{&model_file, &obser_para, &physic_para, &mag_para, &output_file}
	);

	// 查看是否通过强制参数检查
	if (!fp.pass_mandatory()) return 0;

	if (fp.set_opt('g')) components_sign[0] = true;
	if (fp.set_opt('x')) components_sign[1] = true;
	if (fp.set_opt('y')) components_sign[2] = true;
	if (fp.set_opt('z')) components_sign[3] = true;
	if (fp.set_opt('t')) components_sign_mag[0] = true;
	if (fp.set_opt('u')) components_sign_mag[1] = true;
	if (fp.set_opt('v')) components_sign_mag[2] = true;
	if (fp.set_opt('r')) components_sign_mag[3] = true;
	if (fp.set_opt('c')) compact_mode = Packed;

	// check for forward modeling data type
	bool forward_grav = false;
	bool forward_mag = false;

	if (components_sign[0] != false || components_sign[1] != false || 
		components_sign[2] != false || components_sign[3] != false) forward_grav = true;

	if (components_sign_mag[0] != false || components_sign_mag[1] != false || 
		components_sign_mag[2] != false || components_sign_mag[3] != false) forward_mag = true;

	if (forward_grav == false && forward_mag == 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;
	}
	else if (forward_grav == true && forward_mag == true)
	{
		GCTL_ShowWhatError("Can't forward modeling gravitational and magnetic data at the same time.", 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
	int tet_num;
	gctl::array<gctl::vertex3dc> mesh_node;
	gctl::array<gctl::grav_tetrahedron> mesh_tet;
	gctl::array<gctl::gravtet_para> mesh_tet_gp;
	gctl::array<gctl::mag_tetrahedron> mesh_tet_mag;
	gctl::array<gctl::magtet_para> mesh_tet_mp;
	// define physic array
	gctl::array<double> mesh_phys;
	// define result array
	gctl::array<gctl::point3dc> gm_grad;
	gctl::array<gctl::tensor> gm_tensor;

	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.c_str(), "%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.c_str(),"%[^+]+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);
		if (forward_grav) {gctl::read_gmsh_element(infile, mesh_tet, mesh_node, compact_mode); tet_num = mesh_tet.size();}
		if (forward_mag) {gctl::read_gmsh_element(infile, mesh_tet_mag, mesh_node, compact_mode); tet_num = mesh_tet_mag.size();}
		infile.close();
	}
	// proceed as tetgen files
	else
	{
		gctl::read_Tetgen_node(model_file, mesh_node, compact_mode);
		if (forward_grav) {gctl::read_Tetgen_element(model_file, mesh_tet, mesh_node, compact_mode); tet_num = mesh_tet.size();}
		if (forward_mag) {gctl::read_Tetgen_element(model_file, mesh_tet_mag, mesh_node, compact_mode); tet_num = mesh_tet_mag.size();}
	}

	// initiate gravtet_para
	double inc_deg, dec_deg, geoinc_deg, geodec_deg;
	if (forward_grav) gctl::callink_gravity_para(mesh_tet, mesh_tet_gp);
	if (forward_mag)
	{
		gctl::parse_string_to_value(mag_para, '/', true, inc_deg, dec_deg, geoinc_deg, geodec_deg);
		gctl::callink_magnetic_para_earth(mesh_tet_mag, mesh_tet_mp, inc_deg, dec_deg);
	}

	// Firstly try to explain physic_para as float number
	double physic_value;
	if (1 == sscanf(physic_para.c_str(), "%lf", &physic_value))
	{
		mesh_phys.resize(tet_num, 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.c_str(),"%[^+]+d%d", physic_filename, &physic_col))
				strcpy(physic_filename, physic_para.c_str());

			desc.file_name_ = physic_filename;
			gctl::read_text2vector2d(desc, txt_content);

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

			mesh_phys.resize(txt_content.size());
			for (int i = 0; i < mesh_phys.size(); i++)
				mesh_phys.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.c_str());
				strcat(msg, " as a data name in ");
				strcat(msg, model_file.c_str());
				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_phys, physic_para);
		infile.close();
		if (mesh_phys.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 (forward_grav)
	{
		if (components_sign[0])
		{
			gctl::gobser(gm_grad, mesh_tet, obs_points, mesh_phys);
			for (int i = 0; i < obs_points.size(); ++i)
			{
				save_content[i][3] = gm_grad[i].z;
			}

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

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

			if (components_sign[1])
			{
				for (int i = 0; i < obs_points.size(); ++i)
				{
					save_content[i][3] = gm_tensor[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] = gm_tensor[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] = gm_tensor[i].at(2, 2);
				}

				gctl::save_vector2d2text(out_name+"_Vzz.txt", save_content);
			}
		}
	}
	
	if (forward_mag)
	{
		if (components_sign_mag[0])
		{
			gctl::array<double> deltat;
			gctl::magobser(gm_grad, mesh_tet_mag, obs_points, mesh_phys);
			gctl::magnetic_components2deltaT(gm_grad, deltat, geoinc_deg, geodec_deg);

			for (int i = 0; i < obs_points.size(); ++i)
			{
				save_content[i][3] = deltat[i];
			}

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

		if (components_sign_mag[1] || components_sign_mag[2] || components_sign_mag[3])
		{
			gctl::magobser(gm_grad, mesh_tet_mag, obs_points, mesh_phys);

			if (components_sign_mag[1])
			{
				for (int i = 0; i < obs_points.size(); ++i)
				{
					save_content[i][3] = gm_grad[i].x;
				}

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

			if (components_sign_mag[2])
			{
				for (int i = 0; i < obs_points.size(); ++i)
				{
					save_content[i][3] = gm_grad[i].y;
				}

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

			if (components_sign_mag[3])
			{
				for (int i = 0; i < obs_points.size(); ++i)
				{
					save_content[i][3] = gm_grad[i].z;
				}

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