381 lines
14 KiB
C++
381 lines
14 KiB
C++
/********************************************************
|
|
* ██████╗ ██████╗████████╗██╗
|
|
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
|
* ██║ ███╗██║ ██║ ██║
|
|
* ██║ ██║██║ ██║ ██║
|
|
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
|
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
|
* 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::grid_points_2d(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);
|
|
} |