153 lines
7.7 KiB
C++
153 lines
7.7 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.
|
||
******************************************************/
|
||
|
||
#ifndef _GCTL_GRAV_KERNEL_TRICONE_H
|
||
#define _GCTL_GRAV_KERNEL_TRICONE_H
|
||
|
||
#include "gm_data.h"
|
||
|
||
namespace gctl
|
||
{
|
||
struct gravcone_para
|
||
{
|
||
tensor F[4]; ///< tensor products of four facets
|
||
tensor E[12]; ///< tensor products of six edges
|
||
double edglen[12]; ///< edge lengths of six edges
|
||
};
|
||
|
||
typedef type_tricone<gravcone_para> gravcone; ///< 带gravcone_para属性的三角锥结构体
|
||
|
||
/**
|
||
* @brief 基于球面三角剖分的界面反演所使用的三棱锥对象 参考GJI的文章
|
||
*/
|
||
struct gravcone_para_gji : public gravcone_para
|
||
{
|
||
point3dc nface[4];
|
||
point3dc nedge[12];
|
||
//注意三棱锥的侧面与侧面的外法线矢量相对于半径的偏导数为0 所以nf_rj为顶面三角形相对于三个顶点半径的偏导数
|
||
//三棱锥顶面的偏导数矢量
|
||
point3dc nf_rj[3]; //顶面的单位矢量相对于一个标量（顶点半径）的偏导数
|
||
point3dc nj2_rj[3]; //顶面三角形中与顶点相连的逆时针边及其侧面相对于半径的偏导数
|
||
point3dc n2j_rj[3]; //侧面上的偏导数 侧面外法线矢量的偏导数为0
|
||
point3dc n3j_rj[3];
|
||
point3dc nj3_rj[3];
|
||
point3dc n23_rj[3];
|
||
};
|
||
|
||
typedef type_tricone<gravcone_para_gji> gravcone_gji; ///< 带gravcone_para_gji属性的三角锥结构体
|
||
|
||
/**
|
||
* @brief calculate the gravity parameters of given element type.
|
||
*
|
||
* @param in_tri Input elements
|
||
* @param out_para The output parameter array
|
||
*/
|
||
void callink_gravity_para(array<gravcone> &in_cone, array<gravcone_para> &out_para);
|
||
|
||
/**
|
||
* @brief calculate the gravity parameters of given element type.
|
||
*
|
||
* @param in_tri Input elements
|
||
* @param out_para The output parameter array
|
||
*/
|
||
void callink_gravity_para(array<gravcone_gji> &in_cone, array<gravcone_para_gji> &out_para);
|
||
|
||
/**
|
||
* @brief 计算球坐标下的三棱锥单元体正演核矩阵
|
||
*
|
||
* @note 三棱锥类型仅用于球坐标下的重力数据正演
|
||
*
|
||
* @param out_kernel 输出的核矩阵
|
||
* @param[in] ele 三棱锥单元体数组
|
||
* @param[in] ops 球坐标下的观测点数组
|
||
* @param[in] comp_id 待计算的重力场分量 默认为重力值
|
||
* @param[in] show_progress 是否在终端显示进度条 默认为真
|
||
*/
|
||
void gkernel(matrix<double> &out_kernel, const array<gravcone> &ele, const array<point3ds> &obsp,
|
||
gravitational_field_type_e comp_id = Vz, verbose_type_e verbose = FullMsg);
|
||
|
||
/**
|
||
* @brief 计算球坐标下的三棱锥单元体正演重力数据
|
||
*
|
||
* @note 三棱锥类型仅用于球坐标下的重力数据正演
|
||
*
|
||
* @param out_obs 输出的重力数据，重力位的单位为m^2/s^2，重力值单位为m/s^2，重力张量单位为1/s^2
|
||
* @param[in] ele 三棱锥单元体数组，尺度单位为m
|
||
* @param[in] obsp 球坐标下的观测点数组，尺度单位为m
|
||
* @param[in] rho 单元体密度数组，单位为kg/m^3
|
||
* @param[in] comp_id 待计算的重力场分量 默认为重力值
|
||
* @param[in] verbose 返回信息层级
|
||
*/
|
||
void gobser(array<double> &out_obs, const array<gravcone> &ele, const array<point3ds> &obsp,
|
||
const array<double> &rho, gravitational_field_type_e comp_id = Vz, verbose_type_e verbose = FullMsg);
|
||
|
||
/**
|
||
* @brief 计算球坐标下的三棱锥单元体正演重力数据。注意这里计算的是观测点位置x，y和z方向的重力值，而不是球坐标方向的重力值
|
||
*
|
||
* @param out_obs 输出的矢量重力数据，重力位的单位为m^2/s^2，重力值单位为m/s^2，重力张量单位为1/s^2
|
||
* @param[in] ele 三棱锥单元体数组，尺度单位为m
|
||
* @param[in] obsp 球坐标下的观测点数组，尺度单位为m
|
||
* @param[in] rho 单元体密度数组，单位为kg/m^3
|
||
* @param[in] verbose 返回信息层级
|
||
*/
|
||
void gobser(array<point3dc> &out_obs, const array<gravcone> &ele, const array<point3dc> &obsp,
|
||
const array<double> &rho, verbose_type_e verbose = FullMsg);
|
||
|
||
/**
|
||
* @brief 计算球坐标下的三棱锥单元体正演重力数据
|
||
*
|
||
* @note 三棱锥类型仅用于球坐标下的重力数据正演
|
||
*
|
||
* @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<gravcone_gji> &ele, const array<point3ds> &obsp,
|
||
const array<double> &rho, gravitational_field_type_e comp_id = Vz, verbose_type_e verbose = FullMsg);
|
||
|
||
/**
|
||
* @brief 计算球坐标三角剖分下重力观测点位相对于三棱锥单元体半径的导数（用于界面起伏反演 参考GJI文章）
|
||
*
|
||
* @param out_mGrad 输出的模型偏导数
|
||
* @param[in] vert_neighList 三角剖分中每个顶点周围的三棱锥单元体列表
|
||
* @param[in] obsp 观测点数组
|
||
* @param[in] rho 三棱锥密度数组
|
||
* @param[in] verts_ptr 三棱锥顶点数组的指针 这个变量主要用于计算观测点与顶点之间的夹角 从而加速导数的计算 不是必须的
|
||
* @param[in] ang_limit 上一个参数中观测点与顶点之间允许的夹角的最大值 对于某个观测点，超过这个角度的顶点不计算梯度 默认为-1 此时不计算夹角
|
||
* @param[in] comp_id 重力观测值的分量类型 默认为重力值
|
||
* @param[in] verbose 返回信息层级
|
||
*/
|
||
void gobser_model_gradient(array<double> &out_mGrad, const array<std::vector<gravcone_gji*> > &vert_neighList,
|
||
const array<double> &obs_diff, const array<point3ds> &obsp, const array<double> &rho,
|
||
const array<vertex3dc> *verts_ptr = nullptr, double ang_limit = -1.0,
|
||
gravitational_field_type_e comp_id = Vz, verbose_type_e verbose = FullMsg);
|
||
}
|
||
|
||
#endif // _GCTL_GRAV_KERNEL_TRICONE_H
|