tmp update

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)

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 <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"
+
+using namespace gctl;
+
+int main(int argc, char const *argv[]) try
+{
+	// set up observation points
+	array<point2dc> 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<vertex2dc> 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<double> rho(1, 1.0);
+	array<triangle2d> 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<double> 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);
+}

lib/potential/gkernel_triangle2d.cpp

@@ -44,6 +44,10 @@ void gkernel_triangle2d_vz(gctl::matrix<double> &out_kernel, const gctl::array<g
 	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose);
 void gkernel_triangle2d_vx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
 	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose);
+void gkernel_triangle2d_vxx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose);
+void gkernel_triangle2d_vxz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose);
 void gkernel_triangle2d_vzx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
 	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose);
 void gkernel_triangle2d_vzz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
@@ -70,6 +74,12 @@ void gctl::gkernel(matrix<double> &out_kernel, const array<triangle2d> &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<double> &out_kernel, const array<triangle2d> &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<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
+
+void gkernel_triangle2d_vz2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
+void gkernel_triangle2d_vx2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
+void gkernel_triangle2d_vxx2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
+void gkernel_triangle2d_vxz2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
+void gkernel_triangle2d_vzx2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
+void gkernel_triangle2d_vzz2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
+
+void gctl::gkernel(matrix<double> &out_kernel, const array<triangle2d> &ele, 
+	const array<point2dp> &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<double> &out_obs, const gctl::array<gctl::
 	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
 void gobser_triangle2d_vx(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
 	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
+void gobser_triangle2d_vxx(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
+void gobser_triangle2d_vxz(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
 void gobser_triangle2d_vzx(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
 	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
 void gobser_triangle2d_vzz(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
@@ -126,6 +187,12 @@ void gctl::gobser(array<double> &out_obs, const array<triangle2d> &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<double> &out_obs, const array<triangle2d> &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<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
+
+void gobser_triangle2d_vz2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
+void gobser_triangle2d_vx2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
+void gobser_triangle2d_vxx2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
+void gobser_triangle2d_vxz2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
+void gobser_triangle2d_vzx2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
+void gobser_triangle2d_vzz2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
+
+void gctl::gobser(array<double> &out_obs, const array<triangle2d> &ele, const array<point2dp> &obsp, 
+	const array<double> &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<double> &out_kernel, const gctl::array<g
 	return;
 }
 
+void gkernel_triangle2d_vxx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dc> &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<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dc> &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<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
 	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose)
 {
@@ -263,6 +427,190 @@ void gkernel_triangle2d_vzz(gctl::matrix<double> &out_kernel, const gctl::array<
 	return;
 }
 
+void gkernel_triangle2d_vz2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &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<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &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<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &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<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &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<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &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<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &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<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
 	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
 {
@@ -311,6 +659,54 @@ void gobser_triangle2d_vx(gctl::array<double> &out_obs, const gctl::array<gctl::
 	return;
 }
 
+void gobser_triangle2d_vxx(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &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<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &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<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
 	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
 {
@@ -359,6 +755,186 @@ void gobser_triangle2d_vzz(gctl::array<double> &out_obs, const gctl::array<gctl:
 	return;
 }
 
+void gobser_triangle2d_vz2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &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<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &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<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &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<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &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<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &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<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
+	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &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));
 

lib/potential/gkernel_triangle2d.h

@@ -56,6 +56,31 @@ namespace gctl
 	 */
 	void gobser(array<double> &out_obs, const array<triangle2d> &ele, const array<point2dc> &obsp, 
 		const array<double> &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<double> &out_kernel, const array<triangle2d> &ele, const array<point2dp> &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<double> &out_obs, const array<triangle2d> &ele, const array<point2dp> &obsp, 
+		const array<double> &rho, gravitational_field_type_e comp_id = Vz, verbose_type_e verbose = FullMsg);
 }
 
 #endif // _GCTL_GRAV_KERNEL_TRIANGLE2D_H

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"