/********************************************************
* ██████╗ ██████╗████████╗██╗
* ██╔════╝ ██╔════╝╚══██╔══╝██║
* ██║ ███╗██║ ██║ ██║
* ██║ ██║██║ ██║ ██║
* ╚██████╔╝╚██████╗ ██║ ███████╗
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
* Geophysical Computational Tools & Library (GCTL)
*
* Copyright (c) 2023 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.
******************************************************/
#ifndef _GCTL_GEOMETRY3D_H
#define _GCTL_GEOMETRY3D_H
// library's head file
#include "../core/array.h"
#include "../core/matrix.h"
#include "../maths.h"
#include "tensor.h"
#include "node.h"
#include "edge.h"
#include "tetrahedron.h"
#include "triangle.h"
#include "block.h"
#include "tri_cone.h"
#include "sphere.h"
#include "tesseroid.h"
#include "prism.h"
namespace gctl
{
namespace geometry3d
{
/**
* @brief 两个矢量以坐标原点为顶点的夹角（弧度）
*
* @param[in] a 三维空间内的一个向量或实点a。
* @param[in] b 三维空间内的一个向量或实点b。
*
* @return 两个矢量之间的夹角。
*/
double angle(const point3dc &a, const point3dc &b);
/**
* @brief 计算三角形的面积
*
* @param a 三维空间内的一个实点a
* @param b 三维空间内的一个实点b
* @param c 三维空间内的一个实点c
* @return double 面积
*/
double triangle_area(const point3dc &a, const point3dc &b, const point3dc &c);
/**
* @brief 计算四面体的体积
*
* @param a 三维空间内的一个实点a
* @param b 三维空间内的一个实点b
* @param c 三维空间内的一个实点c
* @param d 三维空间内的一个实点d
* @return double 体积
*/
double tetrahedron_volume(const point3dc &a, const point3dc &b, const point3dc &c, const point3dc &d);
/**
* @brief 点到线段的距离
*
* @param[in] line_start 线段的起点
* @param[in] line_end 线段的终点
* @param[in] dot 线段外的一个点
*
* @return 点到线段的距离
*/
double dot2line(const point3dc &line_start, const point3dc &line_end, const point3dc &dot);
/**
* @brief 点在直线上的投影
*
* @param[in] line_start 线段的起点
* @param[in] line_end 线段的终点
* @param[in] dot 线段外的一个点
*
* @return 点在直线上的投影
*/
point3dc dot_on_line(const point3dc &line_start, const point3dc &line_end, const point3dc &dot);
/**
* @brief 点到平面的距离。
*
* @param[in] c 平面上的一个点
* @param[in] n 平面的单位外法线矢量
* @param[in] d 平面外的一点
*
* @return 平面外的点到平面的距离
*/
double dot2plane(const point3dc &c, const point3dc &n, const point3dc &d);
/**
* @brief Calculate the coefficients of a 3D plane using coordinates of three points on the plane.
* The plane's equation is given as Ax + By + Cz + D = 0
*
* @param x1 x coordinate of the first point
* @param x2 x coordinate of the second point
* @param x3 x coordinate of the third point
* @param y1 y coordinate of the first point
* @param y2 y coordinate of the second point
* @param y3 y coordinate of the third point
* @param z1 z coordinate of the first point
* @param z2 z coordinate of the second point
* @param z3 z coordinate of the third point
* @param A coefficient A
* @param B coefficient B
* @param C coefficient C
* @param D coefficient D
*/
void get_plane_coeff(double x1, double x2, double x3, double y1, double y2, double y3,
double z1, double z2, double z3, double &A, double &B, double &C, double &D);
/**
* @brief Calculate the coefficients of a 3D plane defined by two points on a sphercial surface and the origin.
* The plane's equation is given as Ax + By + Cz + D = 0
*
* @param p1 point 1 on the spherical surface.
* @param p2 point 2 on the spherical surface.
* @param A coefficient A
* @param B coefficient B
* @param C coefficient C
* @param D coefficient D
*/
void get_plane_coeff(const point3ds &p1, const point3ds &p2, double &A, double &B, double &C, double &D);
/**
* @brief 计算直线与平面的交点。
*
* @param[in] face_p 平面上的一个点。
* @param[in] face_nor 平面的单位外法线矢量。
* @param[in] line_p 直线上的一个点。
* @param[in] line_nor 直线的单位方向矢量。
*
* @return 返回的矢量。
*/
point3dc line_on_plane(const point3dc &face_p, const point3dc &face_nor,
const point3dc &line_p, const point3dc &line_nor);
/**
* @brief 计算两个矢量之间的椭圆弧上的任意一点
*
* 三维空间中任意两个以原点为起点矢量可视为一个三维椭圆弧上的两个点，给定这两个点与需要计算的第三点矢量
* 与起点矢量的夹角可以获得计算点坐标。
*
* @param[in] v1 起点矢量
* @param[in] v2 终点矢量
* @param[in] phi 计算点矢量与起点矢量之间的夹角(弧度)
*
* @return 返回计算点坐标
*/
point3dc track_ellipse(const point3dc &v1, const point3dc &v2, double phi);
/**
* @brief 点在平面上的投影
*
* @param[in] face_dot 平面上一点
* @param[in] face_nor 平面的法向量
* @param[in] dot 平面外一点
*
* @return 平面外一点在平面上的投影点
*/
point3dc dot_on_plane(const point3dc &face_dot, const point3dc &face_nor, const point3dc &dot);
/**
* @brief 距离平方反比形式的三角形内插值函数
*
* @param[in] p 待插值点坐标
* @param[in] p1 三角形顶点坐标1
* @param[in] p2 三角形顶点坐标2
* @param[in] p3 三角形顶点坐标3
* @param[in] d1 顶点1的值
* @param[in] d2 顶点2的值
* @param[in] d3 顶点3的值
*
* @return 插值点的值
*/
double tri_interp_dis(const point3dc &p, const point3dc &p1, const point3dc &p2,
const point3dc &p3, double d1, double d2, double d3);
/**
* @brief 球心角平方反比形式的三角形内插值函数
*
* @param[in] p 待插值点坐标
* @param[in] p1 三角形顶点坐标1
* @param[in] p2 三角形顶点坐标2
* @param[in] p3 三角形顶点坐标3
* @param[in] d1 顶点1的值
* @param[in] d2 顶点2的值
* @param[in] d3 顶点3的值
*
* @return 插值点的值
*/
double tri_interp_ang(const point3dc &p, const point3dc &p1, const point3dc &p2,
const point3dc &p3, double d1, double d2, double d3);
/**
* @brief 细分三角形面积比形式的三角形内插值函数
*
* @param[in] p 待插值点坐标
* @param[in] p1 三角形顶点坐标1
* @param[in] p2 三角形顶点坐标2
* @param[in] p3 三角形顶点坐标3
* @param[in] d1 顶点1的值
* @param[in] d2 顶点2的值
* @param[in] d3 顶点3的值
*
* @return 插值点的值
*/
double tri_interp_area(const point3dc &p, const point3dc &p1, const point3dc &p2,
const point3dc &p3, double d1, double d2, double d3);
/**
* @brief 距离平方反比形式的三角形内插值函数
*
* @param[in] p 待插值点坐标
* @param[in] p1 三角形顶点坐标1
* @param[in] p2 三角形顶点坐标2
* @param[in] p3 三角形顶点坐标3
* @param[in] p4 三角形顶点坐标4
* @param[in] d1 顶点1的值
* @param[in] d2 顶点2的值
* @param[in] d3 顶点3的值
* @param[in] d4 顶点4的值
*
* @return 插值点的值
*/
double tet_interp_dis(const point3dc &p, const point3dc &p1, const point3dc &p2,
const point3dc &p3, const point3dc &p4, double d1, double d2, double d3, double d4);
/**
* @brief 球坐标转直角坐标
*
* @param[in] sp_ptr 球坐标点的指针
* @param cp_ptr 直角坐标点的指针
*
* @return 执行是否成功
*/
int sph2car(const point3ds *sp_ptr, point3dc *cp_ptr);
/**
* @brief 直角坐标转球坐标
*
* @param[in] cp_ptr 直角坐标点的指针
* @param sp_ptr 球坐标点的指针
*
* @return 执行是否成功
*/
int car2sph(const point3dc *cp_ptr, point3ds *sp_ptr);
/**
* @brief 椭球坐标转直角坐标
*
* @param lon 经度
* @param lat 维度
* @param r 短半径
* @param R 长半径
* @return point3dc 直角坐标系点
*/
point3dc ellip2car(double lon, double lat, double r, double R);
/**
* @brief 判断射线是否过三角形 在order_triangular_surface中调用
*
* @param[in] origin 射线的起点
* @param[in] direct 射线的方向
* @param[in] tri_facet 待检测的三角形
* @param[in] active_edge 长度为3的数组。表示对应三角形的边是否为激活的状态（激活表示如果射线与此边相交也算与三角形相交）
* @param[in] cutoff_limit 计算叉乘的零值精度，小于此值即认为是0
* @param[in] cross_loc 返回交点的位置
*
* @return 是否在三角形内
*/
bool crossed_tri_facet(const point3dc &origin, const point3dc &direct, const triangle &tri_facet,
const bool *active_edge, double cutoff_limit = 1e-10, point3dc *cross_loc = nullptr);
/**
* @brief 监测点t1与t2是否在点p1，p2与p3构成的平面的同一侧
*
* @param p1 平面上一点
* @param p2 平面上一点
* @param p3 平面上一点
* @param t1 待检测点
* @param t2 待检测点
* @return true 在同一侧
* @return false 不在同一侧
*/
bool same_side(const point3dc &p1, const point3dc &p2, const point3dc &p3, const point3dc &t1, const point3dc &t2);
/**
* @brief 判断一个点是否在四面体内
*
* 注意点落在四面体的顶点，边或者面上也算是在四面体内
*
* @param test_p 待检测点
* @param p1 第一个点
* @param p2 第二个点
* @param p3 第三个点
* @param p4 第四个点
*
* @return 是否在四面体内
*/
bool node_in_tetrahedron(const point3dc &p1, const point3dc &p2,
const point3dc &p3, const point3dc &p4, const point3dc &test_p);
/**
* @brief 对多面体表面的三角形的顶点进行排序，保证所有三角形顶点都是逆时针排序的
*
* @param poly_tri_ptr Pointer of the triangle array
*
* @return status of the function
*/
void order_triangular_surface(array &poly_tri, double cutoff_limit = 1e-10);
/**
* @brief 获取四面体网络中的公共面及公共面两侧的四面体指针
*
* @param[in] ele 四面体网络的元素数组
* @param out_face 返回的公共三角形列表数组
* @param out_neigh_ptr 返回的公共三角形两侧邻接四面体的指针数组的指针，
* 若指针不为空则返回邻接四面体指针数组。否则不计算此数组。
*/
void get_tetra_mesh_common_faces(const array &ele, array &out_face,
array *out_neigh_ptr = nullptr);
/**
* @brief 计算四面体网络的顶点数量
*
*/
size_t sort_node_number(const array &ele);
/**
* @brief 计算四面体网络中与顶点相连的四面体的指针列表
*
* @param[in] ele 四面体数组
* @param neigh_list 返回的四面体指针列表
* @param node_nump 四面体数组内顶点个数的指针
*/
void sort_node_neighbor(const array &ele, std::vector > &neigh_list, size_t *node_nump = nullptr);
/**
* @brief 计算四面体网络内的各个四面体之间的相邻关系
*
* @param ele 四面体数组
* @param node_nump 四面体数组内顶点个数的指针
*/
void sort_tet_neighbor(array &ele, size_t *node_nump = nullptr);
/**
* @brief 使用平面切割三角网络
*
* @note 平面法线正方向的三角网络将被保留
*
* @param[in] in_eles 输入的三维三角形数组（一般为成片的三角网）
* @param[in] nor 平面的法线方向
* @param[in] surf 平面上的一点
* @param out_nodes 切割后的新网络顶点集
* @param out_eles 切割后的新网络三角形元素集
*/
void cut_triangular_mesh(const array &in_eles, const point3dc &nor, const point3dc &surf,
array &out_nodes, array &out_eles);
}
}
#endif // _GCTL_GEOMETRY3D_H