diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt
index 17469d0..14588d5 100644
--- a/example/CMakeLists.txt
+++ b/example/CMakeLists.txt
@@ -12,12 +12,13 @@ macro(add_example name switch)
endif()
endmacro()
+add_example(gobser_tri2d_ex ON)
add_example(mobser_dipole_ex OFF)
add_example(mobser_block_ex OFF)
add_example(mobser_block_gradient_ex OFF)
add_example(mobser_tri_ex OFF)
add_example(mobser_tri_sph_ex OFF)
-add_example(mobser_tricone_ex ON)
+add_example(mobser_tricone_ex OFF)
add_example(mobser_tetra_ex OFF)
add_example(mobser_tetra_ex2 OFF)
add_example(mobser_tetra_sph_ex OFF)
diff --git a/example/gobser_tri2d_ex.cpp b/example/gobser_tri2d_ex.cpp
new file mode 100644
index 0000000..2916c19
--- /dev/null
+++ b/example/gobser_tri2d_ex.cpp
@@ -0,0 +1,74 @@
+/********************************************************
+ * ██████╗ ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║ ███╗██║ ██║ ██║
+ * ██║ ██║██║ ██║ ██║
+ * ╚██████╔╝╚██████╗ ██║ ███████╗
+ * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
+ * 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 .
+ *
+ * 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"
+
+using namespace gctl;
+
+int main(int argc, char const *argv[]) try
+{
+ // set up observation points
+ array obs_loc;
+ grid_points_1d(obs_loc, 0.0, 200.0, 2.0, 0.0);
+ int obs_num = obs_loc.size();
+
+ // set nodes locations
+ array node_vert(3);
+ node_vert[0].set(point2dc(90.0, -30.0));
+ node_vert[1].set(point2dc(100.0, -60.0));
+ node_vert[2].set(point2dc(110.0, -10.0));
+
+ // set triangular elements
+ array rho(1, 1.0);
+ array tri(1);
+ // anti-clock wise
+ tri[0].set(node_vert[0], node_vert[1], node_vert[2]); // set triangle's vertex from existing vertex
+
+ // forward modeling
+ array g(obs_num), gx(obs_num), gzx(obs_num), gzz(obs_num);
+ gobser(g, tri, obs_loc, rho, gctl::Vz);
+ gobser(gx, tri, obs_loc, rho, gctl::Vx);
+ gobser(gzx, tri, obs_loc, rho, gctl::Tzx);
+ gobser(gzz, tri, obs_loc, rho, gctl::Tzz);
+
+ geodsv_io fio;
+ fio.init_table(obs_num, 6);
+ fio.set_column_names({"x", "y", "g", "gx", "gzx", "gzz"});
+ fio.fill_column_point2dc(obs_loc, "x", "y");
+ fio.fill_column(g, "g");
+ fio.fill_column(gx, "gx");
+ fio.fill_column(gzx, "gzx");
+ fio.fill_column(gzz, "gzz");
+ fio.save_text("gobser_tri2d_ex");
+ return 0;
+}
+catch (std::exception &e)
+{
+ GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+}
\ No newline at end of file
diff --git a/lib/potential/gkernel_triangle2d.cpp b/lib/potential/gkernel_triangle2d.cpp
index fca6b8e..b80c3a1 100644
--- a/lib/potential/gkernel_triangle2d.cpp
+++ b/lib/potential/gkernel_triangle2d.cpp
@@ -44,6 +44,10 @@ void gkernel_triangle2d_vz(gctl::matrix &out_kernel, const gctl::array &obsp, gctl::verbose_type_e verbose);
void gkernel_triangle2d_vx(gctl::matrix &out_kernel, const gctl::array &ele,
const gctl::array &obsp, gctl::verbose_type_e verbose);
+void gkernel_triangle2d_vxx(gctl::matrix &out_kernel, const gctl::array &ele,
+ const gctl::array &obsp, gctl::verbose_type_e verbose);
+void gkernel_triangle2d_vxz(gctl::matrix &out_kernel, const gctl::array &ele,
+ const gctl::array &obsp, gctl::verbose_type_e verbose);
void gkernel_triangle2d_vzx(gctl::matrix &out_kernel, const gctl::array &ele,
const gctl::array &obsp, gctl::verbose_type_e verbose);
void gkernel_triangle2d_vzz(gctl::matrix &out_kernel, const gctl::array &ele,
@@ -70,6 +74,12 @@ void gctl::gkernel(matrix &out_kernel, const array &ele,
case Vx:
triangle_kernel = gkernel_triangle2d_vx;
break;
+ case Txx:
+ triangle_kernel = gkernel_triangle2d_vxx;
+ break;
+ case Txz:
+ triangle_kernel = gkernel_triangle2d_vxz;
+ break;
case Tzx:
triangle_kernel = gkernel_triangle2d_vzx;
break;
@@ -77,12 +87,59 @@ void gctl::gkernel(matrix &out_kernel, const array &ele,
triangle_kernel = gkernel_triangle2d_vzz;
break;
default:
- triangle_kernel = gkernel_triangle2d_vz;
+ throw std::invalid_argument("Invalid gravitational field type!");
break;
}
return triangle_kernel(out_kernel, ele, obsp, verbose);
}
+typedef void (*gkernel_triangle2d_ptr2)(gctl::matrix &out_kernel, const gctl::array &ele,
+ const gctl::array &obsp, gctl::verbose_type_e verbose);
+
+void gkernel_triangle2d_vz2(gctl::matrix &out_kernel, const gctl::array &ele,
+ const gctl::array &obsp, gctl::verbose_type_e verbose);
+void gkernel_triangle2d_vx2(gctl::matrix &out_kernel, const gctl::array &ele,
+ const gctl::array &obsp, gctl::verbose_type_e verbose);
+void gkernel_triangle2d_vxx2(gctl::matrix &out_kernel, const gctl::array &ele,
+ const gctl::array &obsp, gctl::verbose_type_e verbose);
+void gkernel_triangle2d_vxz2(gctl::matrix &out_kernel, const gctl::array &ele,
+ const gctl::array &obsp, gctl::verbose_type_e verbose);
+void gkernel_triangle2d_vzx2(gctl::matrix &out_kernel, const gctl::array &ele,
+ const gctl::array &obsp, gctl::verbose_type_e verbose);
+void gkernel_triangle2d_vzz2(gctl::matrix &out_kernel, const gctl::array &ele,
+ const gctl::array &obsp, gctl::verbose_type_e verbose);
+
+void gctl::gkernel(matrix &out_kernel, const array &ele,
+ const array &obsp, gravitational_field_type_e comp_id, verbose_type_e verbose)
+{
+ gkernel_triangle2d_ptr2 triangle_kernel2;
+ switch (comp_id)
+ {
+ case Vz:
+ triangle_kernel2 = gkernel_triangle2d_vz2;
+ break;
+ case Vx:
+ triangle_kernel2 = gkernel_triangle2d_vx2;
+ break;
+ case Txx:
+ triangle_kernel2 = gkernel_triangle2d_vxx2;
+ break;
+ case Txz:
+ triangle_kernel2 = gkernel_triangle2d_vxz2;
+ break;
+ case Tzx:
+ triangle_kernel2 = gkernel_triangle2d_vzx2;
+ break;
+ case Tzz:
+ triangle_kernel2 = gkernel_triangle2d_vzz2;
+ break;
+ default:
+ throw std::invalid_argument("Invalid gravitational field type!");
+ break;
+ }
+ return triangle_kernel2(out_kernel, ele, obsp, verbose);
+}
+
/**
* @brief callback interface of the gravitational observations of 2D triangular elements
*
@@ -99,6 +156,10 @@ void gobser_triangle2d_vz(gctl::array &out_obs, const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose);
void gobser_triangle2d_vx(gctl::array &out_obs, const gctl::array &ele,
const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose);
+void gobser_triangle2d_vxx(gctl::array &out_obs, const gctl::array &ele,
+ const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose);
+void gobser_triangle2d_vxz(gctl::array &out_obs, const gctl::array &ele,
+ const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose);
void gobser_triangle2d_vzx(gctl::array &out_obs, const gctl::array &ele,
const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose);
void gobser_triangle2d_vzz(gctl::array &out_obs, const gctl::array &ele,
@@ -126,6 +187,12 @@ void gctl::gobser(array &out_obs, const array &ele, const ar
case Vx:
triangle_obser = gobser_triangle2d_vx;
break;
+ case Txx:
+ triangle_obser = gobser_triangle2d_vxx;
+ break;
+ case Txz:
+ triangle_obser = gobser_triangle2d_vxz;
+ break;
case Tzx:
triangle_obser = gobser_triangle2d_vzx;
break;
@@ -133,14 +200,63 @@ void gctl::gobser(array &out_obs, const array &ele, const ar
triangle_obser = gobser_triangle2d_vzz;
break;
default:
- triangle_obser = gobser_triangle2d_vz;
+ throw std::invalid_argument("Invalid gravitational field type!");
break;
}
return triangle_obser(out_obs, ele, obsp, rho, verbose);
}
+typedef void (*gobser_triangle2d_ptr2)(gctl::array &out_obs, const gctl::array &ele,
+ const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose);
+
+void gobser_triangle2d_vz2(gctl::array &out_obs, const gctl::array &ele,
+ const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose);
+void gobser_triangle2d_vx2(gctl::array &out_obs, const gctl::array &ele,
+ const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose);
+void gobser_triangle2d_vxx2(gctl::array &out_obs, const gctl::array &ele,
+ const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose);
+void gobser_triangle2d_vxz2(gctl::array &out_obs, const gctl::array &ele,
+ const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose);
+void gobser_triangle2d_vzx2(gctl::array &out_obs, const gctl::array &ele,
+ const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose);
+void gobser_triangle2d_vzz2(gctl::array &out_obs, const gctl::array &ele,
+ const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose);
+
+void gctl::gobser(array &out_obs, const array &ele, const array &obsp,
+ const array &rho, gravitational_field_type_e comp_id, verbose_type_e verbose)
+{
+ gobser_triangle2d_ptr2 triangle_obser2;
+ switch (comp_id)
+ {
+ case Vz:
+ triangle_obser2 = gobser_triangle2d_vz2;
+ break;
+ case Vx:
+ triangle_obser2 = gobser_triangle2d_vx2;
+ break;
+ case Txx:
+ triangle_obser2 = gobser_triangle2d_vxx2;
+ break;
+ case Txz:
+ triangle_obser2 = gobser_triangle2d_vxz2;
+ break;
+ case Tzx:
+ triangle_obser2 = gobser_triangle2d_vzx2;
+ break;
+ case Tzz:
+ triangle_obser2 = gobser_triangle2d_vzz2;
+ break;
+ default:
+ throw std::invalid_argument("Invalid gravitational field type!");
+ break;
+ }
+ return triangle_obser2(out_obs, ele, obsp, rho, verbose);
+}
+
double gkernel_triangle2d_vz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
double gkernel_triangle2d_vx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
+double gkernel_triangle2d_vxx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
+double gkernel_triangle2d_vxz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
double gkernel_triangle2d_vzx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
double gkernel_triangle2d_vzz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
@@ -215,6 +331,54 @@ void gkernel_triangle2d_vx(gctl::matrix &out_kernel, const gctl::array &out_kernel, const gctl::array &ele,
+ const gctl::array &obsp, gctl::verbose_type_e verbose)
+{
+ int i, j;
+ int o_size = obsp.size();
+ int e_size = ele.size();
+
+ out_kernel.resize(o_size, e_size);
+
+ gctl::progress_bar bar(o_size, "gkernel_vxx");
+ for (i = 0; i < o_size; i++)
+ {
+ if (verbose == gctl::FullMsg) bar.progressed(i);
+ else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
+
+#pragma omp parallel for private (j) schedule(guided)
+ for (j = 0; j < e_size; j++)
+ {
+ out_kernel[i][j] = gkernel_triangle2d_vxx_sig(ele.get(j), obsp.get(i));
+ }
+ }
+ return;
+}
+
+void gkernel_triangle2d_vxz(gctl::matrix &out_kernel, const gctl::array &ele,
+ const gctl::array &obsp, gctl::verbose_type_e verbose)
+{
+ int i, j;
+ int o_size = obsp.size();
+ int e_size = ele.size();
+
+ out_kernel.resize(o_size, e_size);
+
+ gctl::progress_bar bar(o_size, "gkernel_vxz");
+ for (i = 0; i < o_size; i++)
+ {
+ if (verbose == gctl::FullMsg) bar.progressed(i);
+ else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
+
+#pragma omp parallel for private (j) schedule(guided)
+ for (j = 0; j < e_size; j++)
+ {
+ out_kernel[i][j] = gkernel_triangle2d_vxz_sig(ele.get(j), obsp.get(i));
+ }
+ }
+ return;
+}
+
void gkernel_triangle2d_vzx(gctl::matrix &out_kernel, const gctl::array &ele,
const gctl::array &obsp, gctl::verbose_type_e verbose)
{
@@ -263,6 +427,190 @@ void gkernel_triangle2d_vzz(gctl::matrix &out_kernel, const gctl::array<
return;
}
+void gkernel_triangle2d_vz2(gctl::matrix &out_kernel, const gctl::array &ele,
+ const gctl::array &obsp, gctl::verbose_type_e verbose)
+{
+ int i, j;
+ int o_size = obsp.size();
+ int e_size = ele.size();
+
+ gctl::point2dc obsc, obsg;
+ out_kernel.resize(o_size, e_size);
+
+ gctl::progress_bar bar(o_size, "gkernel_vz");
+ for (i = 0; i < o_size; i++)
+ {
+ if (verbose == gctl::FullMsg) bar.progressed(i);
+ else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
+
+ obsc = obsp[i].p2c();
+
+#pragma omp parallel for private (j, obsg) schedule(guided)
+ for (j = 0; j < e_size; j++)
+ {
+ obsg.x = gkernel_triangle2d_vx_sig(ele.get(j), &obsc);
+ obsg.y = gkernel_triangle2d_vz_sig(ele.get(j), &obsc);
+ out_kernel[i][j] = obsg.y*sin(obsp[i].arc) + obsg.x*cos(obsp[i].arc);
+ }
+ }
+}
+
+void gkernel_triangle2d_vx2(gctl::matrix &out_kernel, const gctl::array &ele,
+ const gctl::array &obsp, gctl::verbose_type_e verbose)
+{
+ int i, j;
+ int o_size = obsp.size();
+ int e_size = ele.size();
+
+ gctl::point2dc obsc, obsg;
+ out_kernel.resize(o_size, e_size);
+
+ gctl::progress_bar bar(o_size, "gkernel_vx");
+ for (i = 0; i < o_size; i++)
+ {
+ if (verbose == gctl::FullMsg) bar.progressed(i);
+ else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
+
+ obsc = obsp[i].p2c();
+
+#pragma omp parallel for private (j, obsg) schedule(guided)
+ for (j = 0; j < e_size; j++)
+ {
+ obsg.x = gkernel_triangle2d_vx_sig(ele.get(j), &obsc);
+ obsg.y = gkernel_triangle2d_vz_sig(ele.get(j), &obsc);
+ out_kernel[i][j] = obsg.y*cos(obsp[i].arc) + obsg.x*sin(obsp[i].arc);
+ }
+ }
+}
+
+void gkernel_triangle2d_vxx2(gctl::matrix &out_kernel, const gctl::array &ele,
+ const gctl::array &obsp, gctl::verbose_type_e verbose)
+{
+ int i, j;
+ int o_size = obsp.size();
+ int e_size = ele.size();
+
+ gctl::point2dc obsc;
+ double t[2][2];
+ out_kernel.resize(o_size, e_size);
+
+ gctl::progress_bar bar(o_size, "gkernel_vxx");
+ for (i = 0; i < o_size; i++)
+ {
+ if (verbose == gctl::FullMsg) bar.progressed(i);
+ else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
+
+ obsc = obsp[i].p2c();
+
+#pragma omp parallel for private (j, t) schedule(guided)
+ for (j = 0; j < e_size; j++)
+ {
+ t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
+ t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
+ t[0][0] = -1.0*t[1][1];
+ t[0][1] = t[1][0];
+ out_kernel[i][j] = cos(obsp[i].arc)*(t[0][0]*cos(obsp[i].arc) + t[1][0]*sin(obsp[i].arc))
+ + sin(obsp[i].arc)*(t[0][1]*cos(obsp[i].arc) + t[1][1]*sin(obsp[i].arc));
+ }
+ }
+}
+
+void gkernel_triangle2d_vxz2(gctl::matrix &out_kernel, const gctl::array &ele,
+ const gctl::array &obsp, gctl::verbose_type_e verbose)
+{
+ int i, j;
+ int o_size = obsp.size();
+ int e_size = ele.size();
+
+ gctl::point2dc obsc;
+ double t[2][2];
+ out_kernel.resize(o_size, e_size);
+
+ gctl::progress_bar bar(o_size, "gkernel_vxz");
+ for (i = 0; i < o_size; i++)
+ {
+ if (verbose == gctl::FullMsg) bar.progressed(i);
+ else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
+
+ obsc = obsp[i].p2c();
+
+#pragma omp parallel for private (j, t) schedule(guided)
+ for (j = 0; j < e_size; j++)
+ {
+ t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
+ t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
+ t[0][0] = -1.0*t[1][1];
+ t[0][1] = t[1][0];
+ out_kernel[i][j] = sin(obsp[i].arc)*(t[0][0]*cos(obsp[i].arc) + t[1][0]*sin(obsp[i].arc))
+ + cos(obsp[i].arc)*(t[0][1]*cos(obsp[i].arc) + t[1][1]*sin(obsp[i].arc));
+ }
+ }
+}
+
+void gkernel_triangle2d_vzx2(gctl::matrix &out_kernel, const gctl::array &ele,
+ const gctl::array &obsp, gctl::verbose_type_e verbose)
+{
+ int i, j;
+ int o_size = obsp.size();
+ int e_size = ele.size();
+
+ gctl::point2dc obsc;
+ double t[2][2];
+ out_kernel.resize(o_size, e_size);
+
+ gctl::progress_bar bar(o_size, "gkernel_vzx");
+ for (i = 0; i < o_size; i++)
+ {
+ if (verbose == gctl::FullMsg) bar.progressed(i);
+ else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
+
+ obsc = obsp[i].p2c();
+
+#pragma omp parallel for private (j, t) schedule(guided)
+ for (j = 0; j < e_size; j++)
+ {
+ t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
+ t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
+ t[0][0] = -1.0*t[1][1];
+ t[0][1] = t[1][0];
+ out_kernel[i][j] = cos(obsp[i].arc)*(t[0][0]*sin(obsp[i].arc) + t[1][0]*cos(obsp[i].arc))
+ + sin(obsp[i].arc)*(t[0][1]*sin(obsp[i].arc) + t[1][1]*cos(obsp[i].arc));
+ }
+ }
+}
+
+void gkernel_triangle2d_vzz2(gctl::matrix &out_kernel, const gctl::array &ele,
+ const gctl::array &obsp, gctl::verbose_type_e verbose)
+{
+ int i, j;
+ int o_size = obsp.size();
+ int e_size = ele.size();
+
+ gctl::point2dc obsc;
+ double t[2][2];
+ out_kernel.resize(o_size, e_size);
+
+ gctl::progress_bar bar(o_size, "gkernel_vzz");
+ for (i = 0; i < o_size; i++)
+ {
+ if (verbose == gctl::FullMsg) bar.progressed(i);
+ else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
+
+ obsc = obsp[i].p2c();
+
+#pragma omp parallel for private (j, t) schedule(guided)
+ for (j = 0; j < e_size; j++)
+ {
+ t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
+ t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
+ t[0][0] = -1.0*t[1][1];
+ t[0][1] = t[1][0];
+ out_kernel[i][j] = sin(obsp[i].arc)*(t[0][0]*sin(obsp[i].arc) + t[1][0]*cos(obsp[i].arc))
+ + cos(obsp[i].arc)*(t[0][1]*sin(obsp[i].arc) + t[1][1]*cos(obsp[i].arc));
+ }
+ }
+}
+
void gobser_triangle2d_vz(gctl::array &out_obs, const gctl::array &ele,
const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose)
{
@@ -311,6 +659,54 @@ void gobser_triangle2d_vx(gctl::array &out_obs, const gctl::array &out_obs, const gctl::array &ele,
+ const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose)
+{
+ int i, j;
+ int o_size = obsp.size();
+ int e_size = ele.size();
+
+ out_obs.resize(o_size, 0.0);
+
+ gctl::progress_bar bar(e_size, "gobser_vxx");
+ for (j = 0; j < e_size; j++)
+ {
+ if (verbose == gctl::FullMsg) bar.progressed(j);
+ else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
+
+#pragma omp parallel for private (i) schedule(guided)
+ for (i = 0; i < o_size; i++)
+ {
+ out_obs[i] += gkernel_triangle2d_vxx_sig(ele.get(j), obsp.get(i)) * rho[j];
+ }
+ }
+ return;
+}
+
+void gobser_triangle2d_vxz(gctl::array &out_obs, const gctl::array &ele,
+ const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose)
+{
+ int i, j;
+ int o_size = obsp.size();
+ int e_size = ele.size();
+
+ out_obs.resize(o_size, 0.0);
+
+ gctl::progress_bar bar(e_size, "gobser_vxz");
+ for (j = 0; j < e_size; j++)
+ {
+ if (verbose == gctl::FullMsg) bar.progressed(j);
+ else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
+
+#pragma omp parallel for private (i) schedule(guided)
+ for (i = 0; i < o_size; i++)
+ {
+ out_obs[i] += gkernel_triangle2d_vxz_sig(ele.get(j), obsp.get(i)) * rho[j];
+ }
+ }
+ return;
+}
+
void gobser_triangle2d_vzx(gctl::array &out_obs, const gctl::array &ele,
const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose)
{
@@ -359,6 +755,186 @@ void gobser_triangle2d_vzz(gctl::array &out_obs, const gctl::array &out_obs, const gctl::array &ele,
+ const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose)
+{
+ int i, j;
+ int o_size = obsp.size();
+ int e_size = ele.size();
+
+ gctl::point2dc obsc, obsg;
+ out_obs.resize(o_size, 0.0);
+
+ gctl::progress_bar bar(e_size, "gobser_vz");
+ for (j = 0; j < e_size; j++)
+ {
+ if (verbose == gctl::FullMsg) bar.progressed(j);
+ else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
+
+#pragma omp parallel for private (i, obsg, obsc) schedule(guided)
+ for (i = 0; i < o_size; i++)
+ {
+ obsc = obsp[i].p2c();
+ obsg.x = gkernel_triangle2d_vx_sig(ele.get(j), &obsc);
+ obsg.y = gkernel_triangle2d_vz_sig(ele.get(j), &obsc);
+ out_obs[i] += (obsg.y*sin(obsp[i].arc) + obsg.x*cos(obsp[i].arc))*rho[j];
+ }
+ }
+ return;
+}
+
+void gobser_triangle2d_vx2(gctl::array &out_obs, const gctl::array &ele,
+ const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose)
+{
+ int i, j;
+ int o_size = obsp.size();
+ int e_size = ele.size();
+
+ gctl::point2dc obsc, obsg;
+ out_obs.resize(o_size, 0.0);
+
+ gctl::progress_bar bar(e_size, "gobser_vx");
+ for (j = 0; j < e_size; j++)
+ {
+ if (verbose == gctl::FullMsg) bar.progressed(j);
+ else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
+
+#pragma omp parallel for private (i, obsg, obsc) schedule(guided)
+ for (i = 0; i < o_size; i++)
+ {
+ obsc = obsp[i].p2c();
+ obsg.x = gkernel_triangle2d_vx_sig(ele.get(j), &obsc);
+ obsg.y = gkernel_triangle2d_vz_sig(ele.get(j), &obsc);
+ out_obs[i] += (obsg.y*cos(obsp[i].arc) + obsg.x*sin(obsp[i].arc))*rho[j];
+ }
+ }
+ return;
+}
+
+void gobser_triangle2d_vxx2(gctl::array &out_obs, const gctl::array &ele,
+ const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose)
+{
+ int i, j;
+ int o_size = obsp.size();
+ int e_size = ele.size();
+
+ gctl::point2dc obsc;
+ double t[2][2];
+ out_obs.resize(o_size, 0.0);
+
+ gctl::progress_bar bar(e_size, "gobser_vxx");
+ for (j = 0; j < e_size; j++)
+ {
+ if (verbose == gctl::FullMsg) bar.progressed(j);
+ else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
+
+#pragma omp parallel for private (i, t, obsc) schedule(guided)
+ for (i = 0; i < o_size; i++)
+ {
+ obsc = obsp[i].p2c();
+ t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
+ t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
+ t[0][0] = -1.0*t[1][1];
+ t[0][1] = t[1][0];
+ out_obs[i] += -1.0*(cos(obsp[i].arc)*(t[0][0]*cos(obsp[i].arc) + t[1][0]*sin(obsp[i].arc))
+ + sin(obsp[i].arc)*(t[0][1]*cos(obsp[i].arc) + t[1][1]*sin(obsp[i].arc)))*rho[j];
+ }
+ }
+}
+
+void gobser_triangle2d_vxz2(gctl::array &out_obs, const gctl::array &ele,
+ const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose)
+{
+ int i, j;
+ int o_size = obsp.size();
+ int e_size = ele.size();
+
+ gctl::point2dc obsc;
+ double t[2][2];
+ out_obs.resize(o_size, 0.0);
+
+ gctl::progress_bar bar(e_size, "gobser_vxz");
+ for (j = 0; j < e_size; j++)
+ {
+ if (verbose == gctl::FullMsg) bar.progressed(j);
+ else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
+
+#pragma omp parallel for private (i, t, obsc) schedule(guided)
+ for (i = 0; i < o_size; i++)
+ {
+ obsc = obsp[i].p2c();
+ t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
+ t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
+ t[0][0] = -1.0*t[1][1];
+ t[0][1] = t[1][0];
+ out_obs[i] += (sin(obsp[i].arc)*(t[0][0]*cos(obsp[i].arc) + t[1][0]*sin(obsp[i].arc))
+ + cos(obsp[i].arc)*(t[0][1]*cos(obsp[i].arc) + t[1][1]*sin(obsp[i].arc)))*rho[j];
+ }
+ }
+}
+
+void gobser_triangle2d_vzx2(gctl::array &out_obs, const gctl::array &ele,
+ const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose)
+{
+ int i, j;
+ int o_size = obsp.size();
+ int e_size = ele.size();
+
+ gctl::point2dc obsc;
+ double t[2][2];
+ out_obs.resize(o_size, 0.0);
+
+ gctl::progress_bar bar(e_size, "gobser_vzx");
+ for (j = 0; j < e_size; j++)
+ {
+ if (verbose == gctl::FullMsg) bar.progressed(j);
+ else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
+
+#pragma omp parallel for private (i, t, obsc) schedule(guided)
+ for (i = 0; i < o_size; i++)
+ {
+ obsc = obsp[i].p2c();
+ t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
+ t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
+ t[0][0] = -1.0*t[1][1];
+ t[0][1] = t[1][0];
+ out_obs[i] += (cos(obsp[i].arc)*(t[0][0]*sin(obsp[i].arc) + t[1][0]*cos(obsp[i].arc))
+ + sin(obsp[i].arc)*(t[0][1]*sin(obsp[i].arc) + t[1][1]*cos(obsp[i].arc)))*rho[j];
+ }
+ }
+}
+
+void gobser_triangle2d_vzz2(gctl::array &out_obs, const gctl::array &ele,
+ const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose)
+{
+ int i, j;
+ int o_size = obsp.size();
+ int e_size = ele.size();
+
+ gctl::point2dc obsc;
+ double t[2][2];
+ out_obs.resize(o_size, 0.0);
+
+ gctl::progress_bar bar(e_size, "gobser_vzz");
+ for (j = 0; j < e_size; j++)
+ {
+ if (verbose == gctl::FullMsg) bar.progressed(j);
+ else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
+
+#pragma omp parallel for private (i, t, obsc) schedule(guided)
+ for (i = 0; i < o_size; i++)
+ {
+ obsc = obsp[i].p2c();
+ t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
+ t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
+ t[0][0] = -1.0*t[1][1];
+ t[0][1] = t[1][0];
+ out_obs[i] += -1.0*(sin(obsp[i].arc)*(t[0][0]*sin(obsp[i].arc) + t[1][0]*cos(obsp[i].arc))
+ + cos(obsp[i].arc)*(t[0][1]*sin(obsp[i].arc) + t[1][1]*cos(obsp[i].arc)))*rho[j];
+ }
+ }
+}
+
double gkernel_triangle2d_vz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr)
{
double sum;
@@ -376,13 +952,15 @@ double gkernel_triangle2d_vz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_p
Ba = atan2(ele_ptr->vert[n]->y - op_ptr->y, ele_ptr->vert[n]->x - op_ptr->x);
Bb = atan2(ele_ptr->vert[(n+1)%3]->y - op_ptr->y, ele_ptr->vert[(n+1)%3]->x - op_ptr->x);
+ if (ele_ptr->vert[n]->y - op_ptr->y > 0.0 && ele_ptr->vert[n]->x - op_ptr->x < 0.0) Ba -= 2*GCTL_Pi;
+ if (ele_ptr->vert[(n+1)%3]->y - op_ptr->y > 0.0 && ele_ptr->vert[(n+1)%3]->x - op_ptr->x < 0.0) Bb -= 2*GCTL_Pi;
B = (ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x)*(Ba - Bb);
Ca = pow(ele_ptr->vert[(n+1)%3]->x - op_ptr->x,2) + pow(ele_ptr->vert[(n+1)%3]->y - op_ptr->y,2);
Cb = pow(ele_ptr->vert[n]->x - op_ptr->x,2) + pow(ele_ptr->vert[n]->y - op_ptr->y,2);
C = 0.5*(ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y)*(log(Ca) - log(Cb));
- sum += A*(B+C);
+ sum += A*(B + C);
}
return -2.0e+8*GCTL_G0*sum;
}
@@ -404,17 +982,29 @@ double gkernel_triangle2d_vx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_p
Ba = atan2(ele_ptr->vert[n]->y - op_ptr->y, ele_ptr->vert[n]->x - op_ptr->x);
Bb = atan2(ele_ptr->vert[(n+1)%3]->y - op_ptr->y, ele_ptr->vert[(n+1)%3]->x - op_ptr->x);
+ if (ele_ptr->vert[n]->y - op_ptr->y > 0.0 && ele_ptr->vert[n]->x - op_ptr->x < 0.0) Ba -= 2*GCTL_Pi;
+ if (ele_ptr->vert[(n+1)%3]->y - op_ptr->y > 0.0 && ele_ptr->vert[(n+1)%3]->x - op_ptr->x < 0.0) Bb -= 2*GCTL_Pi;
B = (ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y)*(Ba - Bb);
Ca = pow(ele_ptr->vert[(n+1)%3]->x - op_ptr->x,2) + pow(ele_ptr->vert[(n+1)%3]->y - op_ptr->y,2);
Cb = pow(ele_ptr->vert[n]->x - op_ptr->x,2) + pow(ele_ptr->vert[n]->y - op_ptr->y,2);
- C = 0.5*(ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x)*(log(Ca) - log(Cb));
+ C = 0.5*(ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x)*(log(Cb) - log(Ca));
- sum += A*(B+C);
+ sum += A*(B + C);
}
return -2.0e+8*GCTL_G0*sum;
}
+double gkernel_triangle2d_vxx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr)
+{
+ return -1.0*gkernel_triangle2d_vzz_sig(ele_ptr, op_ptr);
+}
+
+double gkernel_triangle2d_vxz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr)
+{
+ return gkernel_triangle2d_vzx_sig(ele_ptr, op_ptr);
+}
+
double gkernel_triangle2d_vzx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr)
{
double sum;
@@ -435,10 +1025,12 @@ double gkernel_triangle2d_vzx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_
Ax= Aa_x/Ab;
Ba = atan2(ele_ptr->vert[n]->y - op_ptr->y, ele_ptr->vert[n]->x - op_ptr->x);
+ if (ele_ptr->vert[n]->y - op_ptr->y > 0.0 && ele_ptr->vert[n]->x - op_ptr->x < 0.0) Ba -= 2*GCTL_Pi;
Ba_x = (ele_ptr->vert[n]->y - op_ptr->y)/((ele_ptr->vert[n]->x - op_ptr->x)*(ele_ptr->vert[n]->x - op_ptr->x)
+ (ele_ptr->vert[n]->y - op_ptr->y)*(ele_ptr->vert[n]->y - op_ptr->y));
Bb = atan2(ele_ptr->vert[(n+1)%3]->y - op_ptr->y, ele_ptr->vert[(n+1)%3]->x - op_ptr->x);
+ if (ele_ptr->vert[(n+1)%3]->y - op_ptr->y > 0.0 && ele_ptr->vert[(n+1)%3]->x - op_ptr->x < 0.0) Bb -= 2*GCTL_Pi;
Bb_x = (ele_ptr->vert[(n+1)%3]->y - op_ptr->y)/((ele_ptr->vert[(n+1)%3]->x - op_ptr->x)*(ele_ptr->vert[(n+1)%3]->x - op_ptr->x)
+ (ele_ptr->vert[(n+1)%3]->y - op_ptr->y)*(ele_ptr->vert[(n+1)%3]->y - op_ptr->y));
@@ -478,10 +1070,12 @@ double gkernel_triangle2d_vzz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_
Ay= Aa_y/Ab;
Ba = atan2(ele_ptr->vert[n]->y - op_ptr->y, ele_ptr->vert[n]->x - op_ptr->x);
+ if (ele_ptr->vert[n]->y - op_ptr->y > 0.0 && ele_ptr->vert[n]->x - op_ptr->x < 0.0) Ba -= 2*GCTL_Pi;
Ba_y = -1.0*(ele_ptr->vert[n]->x - op_ptr->x)/((ele_ptr->vert[n]->x - op_ptr->x)*(ele_ptr->vert[n]->x - op_ptr->x)
+ (ele_ptr->vert[n]->y - op_ptr->y)*(ele_ptr->vert[n]->y - op_ptr->y));
Bb = atan2(ele_ptr->vert[(n+1)%3]->y - op_ptr->y, ele_ptr->vert[(n+1)%3]->x - op_ptr->x);
+ if (ele_ptr->vert[(n+1)%3]->y - op_ptr->y > 0.0 && ele_ptr->vert[(n+1)%3]->x - op_ptr->x < 0.0) Bb -= 2*GCTL_Pi;
Bb_y = -1.0*(ele_ptr->vert[(n+1)%3]->x - op_ptr->x)/((ele_ptr->vert[(n+1)%3]->x - op_ptr->x)*(ele_ptr->vert[(n+1)%3]->x - op_ptr->x)
+ (ele_ptr->vert[(n+1)%3]->y - op_ptr->y)*(ele_ptr->vert[(n+1)%3]->y - op_ptr->y));
diff --git a/lib/potential/gkernel_triangle2d.h b/lib/potential/gkernel_triangle2d.h
index b688403..08079c9 100644
--- a/lib/potential/gkernel_triangle2d.h
+++ b/lib/potential/gkernel_triangle2d.h
@@ -56,6 +56,31 @@ namespace gctl
*/
void gobser(array &out_obs, const array &ele, const array &obsp,
const array &rho, gravitational_field_type_e comp_id = Vz, verbose_type_e verbose = FullMsg);
+
+ /**
+ * @brief 计算极坐标下的二维三角形单元体正演核矩阵
+ *
+ * @param out_kernel 输出的核矩阵
+ * @param[in] ele 二维三角形单元体数组
+ * @param[in] obsp 二维极坐标下的观测点数组
+ * @param[in] comp_id 待计算的重力场分量 默认为重力值
+ * @param[in] verbose 返回信息层级
+ */
+ void gkernel(matrix &out_kernel, const array &ele, const array &obsp,
+ gravitational_field_type_e comp_id = Vz, verbose_type_e verbose = FullMsg);
+
+ /**
+ * @brief 计算极坐标下的二维三角形单元体重力观测值正演
+ *
+ * @param out_obs 输出的重力观测值
+ * @param[in] ele 二维三角形单元体数组
+ * @param[in] obsp 二维极坐标下的观测点数组
+ * @param[in] rho 三角形单元体密度数组
+ * @param[in] comp_id 待计算的重力场分量 默认为重力值
+ * @param[in] verbose 返回信息层级
+ */
+ void gobser(array &out_obs, const array &ele, const array &obsp,
+ const array &rho, gravitational_field_type_e comp_id = Vz, verbose_type_e verbose = FullMsg);
}
#endif // _GCTL_GRAV_KERNEL_TRIANGLE2D_H
\ No newline at end of file
diff --git a/lib/potential/gm_data.h b/lib/potential/gm_data.h
index c48d71c..ad6ea1f 100644
--- a/lib/potential/gm_data.h
+++ b/lib/potential/gm_data.h
@@ -31,6 +31,7 @@
#include "gctl/core.h"
#include "gctl/utility.h"
#include "gctl/geometry.h"
+#include "gctl/maths.h"
#include "gctl/algorithm.h"
#include "gctl_potential_config.h"