139 lines
4.6 KiB
C++
139 lines
4.6 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.
|
|
******************************************************/
|
|
|
|
#include "gctl/core.h"
|
|
#include "gctl/io.h"
|
|
#include "gctl/potential.h"
|
|
#include "iostream"
|
|
|
|
int main(int argc, char *argv[])
|
|
{
|
|
gctl::array<gctl::vertex3dc> msh_vert;
|
|
gctl::array<gctl::grav_tri_cone> msh_cone;
|
|
gctl::array<gctl::gravcone_para> cone_para;
|
|
|
|
// read element from Gmsh file
|
|
std::ifstream infile;
|
|
gctl::open_infile(infile, "data/gobser_tricone/template", ".msh");
|
|
gctl::read_gmsh_node(infile, msh_vert, gctl::Packed);
|
|
gctl::read_gmsh_element(infile, msh_cone, msh_vert, gctl::Packed, nullptr);
|
|
infile.close();
|
|
|
|
// set origin for the cone
|
|
gctl::vertex3dc ori_vert(gctl::point3dc(0.0, 0.0, 0.0), msh_vert.size());
|
|
// create grav_tri_cone
|
|
for (int i = 0; i < msh_cone.size(); i++)
|
|
{
|
|
msh_cone[i].set_origin(ori_vert);
|
|
}
|
|
|
|
//initialize the tensors
|
|
gctl::callink_gravity_para(msh_cone, cone_para);
|
|
|
|
// read topography data
|
|
gctl::_2d_vector txt_content;
|
|
gctl::array<gctl::point3ds> topo_sph;
|
|
gctl::text_descriptor desc;
|
|
gctl::read_text2vector2d("data/gobser_tricone/topo.txt", txt_content, desc);
|
|
|
|
topo_sph.resize(txt_content.size());
|
|
for (int i = 0; i < txt_content.size(); i++)
|
|
{
|
|
topo_sph[i].rad = 1e+5 + txt_content[i][2];
|
|
topo_sph[i].lon = txt_content[i][0];
|
|
topo_sph[i].lat = txt_content[i][1];
|
|
}
|
|
|
|
// initiate observation points
|
|
double xmin = 30, xmax = 50, dx = 0.25;
|
|
double ymin = 30, ymax = 50, dy = 0.25;
|
|
int xnum = round((xmax - xmin)/dx) + 1;
|
|
int ynum = round((ymax - ymin)/dy) + 1;
|
|
gctl::array<gctl::point3ds> obs_ps(xnum*ynum);
|
|
for (int j = 0; j < ynum; j++)
|
|
{
|
|
for (int i = 0; i < xnum; i++)
|
|
{
|
|
obs_ps[i + j*xnum].rad = 1e+5 + 200.0;
|
|
obs_ps[i + j*xnum].lon = xmin + i*dx;
|
|
obs_ps[i + j*xnum].lat = ymin + j*dy;
|
|
}
|
|
}
|
|
|
|
// calculate base gravity
|
|
gctl::array<double> out_obs;
|
|
gctl::array<double> rho(msh_cone.size(), 1.0);
|
|
gctl::gobser(out_obs, msh_cone, obs_ps, rho, gctl::Vz);
|
|
|
|
// interpolate topography
|
|
int m, n;
|
|
double colat1, colat2, lon1, lon2;
|
|
gctl::point3ds tmp_p;
|
|
gctl::point3dc tmp_c;
|
|
|
|
gctl::array<double> out_topo(msh_vert.size());
|
|
for (int i = 0; i < msh_vert.size(); i++)
|
|
{
|
|
tmp_p = msh_vert[i].c2s();
|
|
m = floor(tmp_p.lon - 28.0)/0.5;
|
|
n = floor(tmp_p.lat - 28.0)/0.5;
|
|
|
|
colat1 = 90.0 - (28.0 + n*0.5);
|
|
colat2 = 90.0 - (28.0 + (n+1)*0.5);
|
|
lon1 = 28.0 + m*0.5;
|
|
lon2 = 28.0 + (m+1)*0.5;
|
|
|
|
tmp_p.rad = gctl::sph_linear_interpolate_deg(colat1, colat2, lon1, lon2,
|
|
90.0-tmp_p.lat, tmp_p.lon,
|
|
topo_sph[m + 49*n].rad,
|
|
topo_sph[m + 1 + 49*n].rad,
|
|
topo_sph[m + 49*(n+1)].rad,
|
|
topo_sph[m + 1 + 49*(n+1)].rad);
|
|
|
|
out_topo[i] = tmp_p.rad - 1e+5;
|
|
|
|
tmp_c = tmp_p.s2c();
|
|
msh_vert[i].set(tmp_c);
|
|
}
|
|
|
|
// recalculate tensors
|
|
gctl::callink_gravity_para(msh_cone, cone_para);
|
|
|
|
gctl::array<double> out_obs2;
|
|
gctl::gobser(out_obs2, msh_cone, obs_ps, rho, gctl::Vz);
|
|
|
|
for (int i = 0; i < out_obs.size(); i++)
|
|
{
|
|
out_obs2[i] -= out_obs[i];
|
|
}
|
|
|
|
for (int i = 0; i < out_obs.size(); i++)
|
|
{
|
|
std::cout << obs_ps[i].lon << " " << obs_ps[i].lat << " " << out_obs2[i] << std::endl;
|
|
}
|
|
return 0;
|
|
} |