523 lines
14 KiB
C++
523 lines
14 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 "fmm_data.h"
|
|
|
|
gctl::fmm_triangle2d::fmm_triangle2d()
|
|
{
|
|
reset();
|
|
}
|
|
|
|
void gctl::fmm_triangle2d::reset()
|
|
{
|
|
entity<fmm_vertex2dc, 3>::reset();
|
|
|
|
initialized = false;
|
|
tag = 0;
|
|
slow_ptr = nullptr;
|
|
area = edge_len[0] = edge_len[1] = edge_len[2] = 0.0;
|
|
}
|
|
|
|
/**
|
|
* @brief 初始化结构体
|
|
*/
|
|
void gctl::fmm_triangle2d::init()
|
|
{
|
|
// calculate edge lengths
|
|
for (int i = 0; i < 3; i++)
|
|
{
|
|
// opposite edge length of vertex j
|
|
edge_len[i] = distance(*vert[(i+1)%3], *vert[(i+2)%3]);
|
|
}
|
|
|
|
// calculate element's area
|
|
double p = 0.5*(edge_len[0] + edge_len[1] + edge_len[2]);
|
|
area = sqrt(p*(p-edge_len[0])*(p-edge_len[1])*(p-edge_len[2]));
|
|
|
|
std::vector<fmm_vertex2dc*>::iterator pos;
|
|
for (int i = 0; i < 3; i++)
|
|
{
|
|
// add current element to each vertex element's neighboring list
|
|
vert[i]->e_neigh.push_back(this);
|
|
vert[i]->ve_order.push_back(i);
|
|
// add the other two vertice to current vertex's node neighboring list
|
|
// if the pointer can not be found. add it to the list.
|
|
if (vert[i]->v_neigh.empty())
|
|
vert[i]->v_neigh.push_back(vert[(i+1)%3]);
|
|
else
|
|
{
|
|
pos = std::find(vert[i]->v_neigh.begin(), vert[i]->v_neigh.end(), vert[(i+1)%3]);
|
|
if (pos == vert[i]->v_neigh.end())
|
|
vert[i]->v_neigh.push_back(vert[(i+1)%3]);
|
|
}
|
|
|
|
pos = std::find(vert[i]->v_neigh.begin(), vert[i]->v_neigh.end(), vert[(i+2)%3]);
|
|
if (pos == vert[i]->v_neigh.end())
|
|
vert[i]->v_neigh.push_back(vert[(i+2)%3]);
|
|
}
|
|
|
|
initialized = true;
|
|
return;
|
|
}
|
|
|
|
/**
|
|
* @brief 初始化结构体
|
|
*
|
|
* 此函数会计算三角形的面积、边长、慢度等参数，同时也会
|
|
* 对顶点的v_neigh、e_neigh和ve_order向量进行赋值。
|
|
*
|
|
* @param index 三角形序号
|
|
* @param vec_p0 顶点1指针
|
|
* @param vec_p1 顶点2指针
|
|
* @param vec_p2 顶点3指针
|
|
*/
|
|
void gctl::fmm_triangle2d::init(int index, fmm_vertex2dc *vec_p0,
|
|
fmm_vertex2dc *vec_p1, fmm_vertex2dc *vec_p2, double *slow_p)
|
|
{
|
|
if (index < 0)
|
|
throw invalid_argument("Invalid index. From gctl::fmm_triangle2d::init()");
|
|
|
|
if (vec_p0 == nullptr || vec_p1 == nullptr || vec_p2 == nullptr || slow_p == nullptr)
|
|
throw domain_error("Invalid pointer. From gctl::fmm_triangle2d::init()");
|
|
|
|
id = index;
|
|
slow_ptr = slow_p;
|
|
vert[0] = vec_p0; vert[1] = vec_p1; vert[2] = vec_p2;
|
|
|
|
init();
|
|
return;
|
|
}
|
|
|
|
gctl::fmm_vertex2dc::fmm_vertex2dc() : vertex<point2dc>::vertex()
|
|
{
|
|
reset();
|
|
}
|
|
|
|
void gctl::fmm_vertex2dc::reset()
|
|
{
|
|
tag = 0;
|
|
jn_id1 = jn_id2 = jn_ele = -1;
|
|
jn_grad = jn_ksi = 0.0;
|
|
time_ptr = nullptr;
|
|
if (!v_neigh.empty()) v_neigh.clear();
|
|
if (!e_neigh.empty()) e_neigh.clear();
|
|
if (!ve_order.empty()) ve_order.clear();
|
|
return;
|
|
}
|
|
|
|
gctl::seis_point2d::seis_point2d() : vertex<point2dc>::vertex()
|
|
{
|
|
reset();
|
|
}
|
|
|
|
void gctl::seis_point2d::reset()
|
|
{
|
|
time = GCTL_BDL_MAX;
|
|
return;
|
|
}
|
|
|
|
gctl::seis_point2d_tri::seis_point2d_tri() : seis_point2d::seis_point2d()
|
|
{
|
|
reset();
|
|
}
|
|
|
|
void gctl::seis_point2d_tri::reset()
|
|
{
|
|
host_ele = nullptr;
|
|
return;
|
|
}
|
|
|
|
/**
|
|
* @brief 定位自由点所在的三角形
|
|
*
|
|
* @param in_tris 网格三角形数组指针
|
|
* @param[in] in_num 三角形数组大小
|
|
*/
|
|
void gctl::seis_point2d_tri::find_host_element(fmm_triangle2d *in_tris, int in_num)
|
|
{
|
|
if (in_tris == nullptr)
|
|
throw domain_error("Invalid pointer. From gctl::seis_point2d_tri::find_host_element()");
|
|
|
|
if (in_num <= 0)
|
|
throw invalid_argument("Invalid mesh size. From gctl::seis_point2d_tri::find_host_element()");
|
|
|
|
host_ele = nullptr;
|
|
|
|
bool in_ele;
|
|
point2dc tmp_p[3];
|
|
point2dc l1, l2;
|
|
for (int i = 0; i < in_num; i++)
|
|
{
|
|
if (!in_tris[i].initialized)
|
|
throw runtime_error("Element initialized. From gctl::seis_point2d_tri::find_host_element()");
|
|
|
|
in_ele = true;
|
|
for (int j = 0; j < 3; j++)
|
|
{
|
|
tmp_p[j] = *in_tris[i].vert[j];
|
|
}
|
|
|
|
if ((tmp_p[1]-tmp_p[0]).x*(tmp_p[2]-tmp_p[0]).y -
|
|
(tmp_p[1]-tmp_p[0]).y*(tmp_p[2]-tmp_p[0]).x > 0)
|
|
{
|
|
for (int j = 0; j < 3; j++)
|
|
{
|
|
l1 = tmp_p[(j+1)%3] - tmp_p[j];
|
|
l2.x = x - tmp_p[j].x;
|
|
l2.y = y - tmp_p[j].y;
|
|
// 点落在线上时也算在三角形内
|
|
if (l1.x*l2.y-l1.y*l2.x < 0)
|
|
{
|
|
in_ele = false;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
for (int j = 0; j < 3; j++)
|
|
{
|
|
l1 = tmp_p[(j+1)%3] - tmp_p[j];
|
|
l2.x = x - tmp_p[j].x;
|
|
l2.y = y - tmp_p[j].y;
|
|
if (l1.x*l2.y-l1.y*l2.x > 0)
|
|
{
|
|
in_ele = false;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
if (in_ele)
|
|
{
|
|
host_ele = &in_tris[i];
|
|
break;
|
|
}
|
|
}
|
|
|
|
if (host_ele == nullptr)
|
|
throw runtime_error("No host element found. From gctl::seis_point2d_tri::find_host_element()");
|
|
return;
|
|
}
|
|
|
|
gctl::fmm_tetrahedron::fmm_tetrahedron()
|
|
{
|
|
reset();
|
|
}
|
|
|
|
void gctl::fmm_tetrahedron::reset()
|
|
{
|
|
entity<fmm_vertex3dc, 4>::reset();
|
|
|
|
initialized = false;
|
|
tag = 0;
|
|
slow_ptr = nullptr;
|
|
volume = 0.0;
|
|
face_area[0] = face_area[1] = face_area[2] = face_area[3] = 0.0;
|
|
edge_len[0] = edge_len[1] = edge_len[2] = 0.0;
|
|
edge_len[3] = edge_len[4] = edge_len[5] = 0.0;
|
|
return;
|
|
}
|
|
|
|
/**
|
|
* @brief 初始化结构体
|
|
*/
|
|
void gctl::fmm_tetrahedron::init()
|
|
{
|
|
// calculate edge lengths
|
|
int count = 0;
|
|
// 01 02 03 12 13 23
|
|
for (int i = 0; i < 3; i++)
|
|
{
|
|
for (int j = i+1; j < 4; j++)
|
|
{
|
|
edge_len[count] = distance(*vert[i], *vert[j]);
|
|
count++;
|
|
}
|
|
}
|
|
|
|
// calculate element's area
|
|
double p;
|
|
p = 0.5*(edge_len[0] + edge_len[1] + edge_len[3]); // 012
|
|
face_area[0] = sqrt(p*(p-edge_len[0])*(p-edge_len[1])*(p-edge_len[3]));
|
|
|
|
p = 0.5*(edge_len[0] + edge_len[2] + edge_len[4]); // 013
|
|
face_area[1] = sqrt(p*(p-edge_len[0])*(p-edge_len[2])*(p-edge_len[4]));
|
|
|
|
p = 0.5*(edge_len[1] + edge_len[2] + edge_len[5]); // 023
|
|
face_area[2] = sqrt(p*(p-edge_len[1])*(p-edge_len[2])*(p-edge_len[5]));
|
|
|
|
p = 0.5*(edge_len[3] + edge_len[4] + edge_len[5]); // 123
|
|
face_area[3] = sqrt(p*(p-edge_len[3])*(p-edge_len[4])*(p-edge_len[5]));
|
|
|
|
// calculate element's volume
|
|
gctl::point3dc va, vb, vc;
|
|
va = *vert[1] - *vert[0];
|
|
vb = *vert[2] - *vert[0];
|
|
vc = *vert[3] - *vert[0];
|
|
volume = GCTL_FABS(dot(cross(va,vb),vc))/6.0;
|
|
|
|
std::vector<fmm_vertex3dc*>::iterator pos;
|
|
for (int i = 0; i < 4; i++)
|
|
{
|
|
// add current element to each vertex element's neighboring list
|
|
vert[i]->e_neigh.push_back(this);
|
|
vert[i]->ve_order.push_back(i);
|
|
// add the other three vertice to current vertex's node neighboring list
|
|
// if the pointer can not be found. add it to the list.
|
|
if (vert[i]->v_neigh.empty())
|
|
vert[i]->v_neigh.push_back(vert[(i+1)%4]);
|
|
else
|
|
{
|
|
pos = std::find(vert[i]->v_neigh.begin(), vert[i]->v_neigh.end(),
|
|
vert[(i+1)%4]);
|
|
if (pos == vert[i]->v_neigh.end())
|
|
vert[i]->v_neigh.push_back(vert[(i+1)%4]);
|
|
}
|
|
|
|
pos = std::find(vert[i]->v_neigh.begin(), vert[i]->v_neigh.end(),
|
|
vert[(i+2)%4]);
|
|
if (pos == vert[i]->v_neigh.end())
|
|
vert[i]->v_neigh.push_back(vert[(i+2)%4]);
|
|
|
|
pos = std::find(vert[i]->v_neigh.begin(), vert[i]->v_neigh.end(),
|
|
vert[(i+3)%4]);
|
|
if (pos == vert[i]->v_neigh.end())
|
|
vert[i]->v_neigh.push_back(vert[(i+3)%4]);
|
|
}
|
|
|
|
initialized = true;
|
|
return;
|
|
}
|
|
|
|
/**
|
|
* @brief 初始化结构体
|
|
*
|
|
* 此函数会计算四面体的体积、三角形面的面积、四面体边长、慢度等参数，同时也会
|
|
* 对顶点的v_neigh、e_neigh和ve_order向量进行赋值。
|
|
*
|
|
* @param index 三角形序号
|
|
* @param vec_p0 顶点1指针
|
|
* @param vec_p1 顶点2指针
|
|
* @param vec_p2 顶点3指针
|
|
* @param vec_p3 顶点4指针
|
|
*/
|
|
void gctl::fmm_tetrahedron::init(int index, fmm_vertex3dc *vec_p0,
|
|
fmm_vertex3dc *vec_p1, fmm_vertex3dc *vec_p2, fmm_vertex3dc *vec_p3)
|
|
{
|
|
if (index < 0)
|
|
throw invalid_argument("Invalid index. From gctl::fmm_tetrahedron::init()");
|
|
|
|
if (vec_p0 == nullptr || vec_p1 == nullptr ||
|
|
vec_p2 == nullptr || vec_p3 == nullptr)
|
|
throw domain_error("Invalid pointer. From gctl::fmm_tetrahedron::init()");
|
|
|
|
id = index;
|
|
vert[0] = vec_p0; vert[1] = vec_p1;
|
|
vert[2] = vec_p2; vert[3] = vec_p3;
|
|
|
|
init();
|
|
return;
|
|
}
|
|
|
|
int gctl::fmm_tetrahedron::vert2edge_id(int v1d, int v2d)
|
|
{
|
|
if (GCTL_MIN(v1d, v2d) == 0)
|
|
return GCTL_MAX(v1d, v2d)-1; // 01->0 02->1 03->2
|
|
else
|
|
return v1d+v2d; // 12->3 13->4 23->5
|
|
}
|
|
|
|
int gctl::fmm_tetrahedron::vert2face_id(int v1d, int v2d, int v3d)
|
|
{
|
|
if (GCTL_MIN(GCTL_MIN(v1d, v2d), v3d) == 1)
|
|
return 3; // 123
|
|
else if (GCTL_MAX(GCTL_MAX(v1d, v2d), v3d) == 2)
|
|
return 0; // 012
|
|
else if (v1d + v2d + v3d == 4)
|
|
return 1; // 013
|
|
else
|
|
return 2; // 023
|
|
}
|
|
|
|
gctl::fmm_vertex3dc::fmm_vertex3dc() : vertex<point3dc>::vertex()
|
|
{
|
|
reset();
|
|
}
|
|
|
|
void gctl::fmm_vertex3dc::reset()
|
|
{
|
|
tag = 0;
|
|
jn_id1 = jn_id2 = jn_id3 = jn_ele = -1;
|
|
jn_grad = jn_ksi = jn_zeta = 0.0;
|
|
time_ptr = nullptr;
|
|
if (!v_neigh.empty()) v_neigh.clear();
|
|
if (!e_neigh.empty()) e_neigh.clear();
|
|
if (!ve_order.empty()) ve_order.clear();
|
|
return;
|
|
}
|
|
|
|
gctl::seis_point3d::seis_point3d() : vertex<point3dc>::vertex()
|
|
{
|
|
reset();
|
|
}
|
|
|
|
void gctl::seis_point3d::reset()
|
|
{
|
|
time = GCTL_BDL_MAX;
|
|
return;
|
|
}
|
|
|
|
gctl::seis_point3d_tet::seis_point3d_tet()
|
|
{
|
|
reset();
|
|
}
|
|
|
|
void gctl::seis_point3d_tet::reset()
|
|
{
|
|
host_ele = nullptr;
|
|
return;
|
|
}
|
|
|
|
/**
|
|
* @brief 定位自由点所在的四面体
|
|
*
|
|
* @param in_tets 网格四面体数组指针
|
|
* @param[in] in_num 三角形数组大小
|
|
*/
|
|
void gctl::seis_point3d_tet::find_host_element(fmm_tetrahedron *in_tets, int in_num)
|
|
{
|
|
if (in_tets == nullptr)
|
|
throw domain_error("Invalid pointer. gctl::seis_point3d_tet::find_host_element()");
|
|
|
|
if (in_num <= 0)
|
|
throw invalid_argument("Invalid mesh size. From gctl::seis_point3d_tet::find_host_element()");
|
|
|
|
host_ele = nullptr;
|
|
|
|
point3dc tmp(x, y, z);
|
|
for (int i = 0; i < in_num; i++)
|
|
{
|
|
if (!in_tets[i].initialized)
|
|
throw runtime_error("Element initialized. From gctl::seis_point3d_tet::find_host_element()");
|
|
|
|
if (geometry3d::node_in_tetrahedron(*in_tets[i].vert[0],
|
|
*in_tets[i].vert[1], *in_tets[i].vert[2], *in_tets[i].vert[3], tmp))
|
|
{
|
|
host_ele = &in_tets[i];
|
|
break;
|
|
}
|
|
}
|
|
|
|
if (host_ele == nullptr)
|
|
throw runtime_error("No host element found for location (" + std::to_string(x) + ", " + std::to_string(y) + ", " + std::to_string(z) + "). From gctl::seis_point3d_tet::find_host_element()");
|
|
return;
|
|
}
|
|
|
|
/**
|
|
* @brief 创建一个二维三角网络的FMM正演网络
|
|
*
|
|
* @param[in] in_node 输入的二维顶点集
|
|
* @param[in] in_ele 输入的二维三角形集
|
|
* @param[in] in_time 输入的顶点顶点走时
|
|
* @param out_node 输出的二维FMM顶点集
|
|
* @param out_ele 输出的二维三角形集
|
|
*/
|
|
void gctl::create_fmm_mesh(const array<vertex2dc> &in_node, const array<triangle2d> &in_ele,
|
|
const array<double> &in_time, const array<double> &in_slow, array<fmm_vertex2dc> &out_node,
|
|
array<fmm_triangle2d> &out_ele)
|
|
{
|
|
if (in_node.empty() || in_ele.empty() || in_time.empty())
|
|
throw runtime_error("Invalid array sizes. From void gctl::create_fmm_mesh(...)");
|
|
|
|
if (in_node.size() != in_time.size())
|
|
throw runtime_error("Incompatible array sizes. From void gctl::create_fmm_mesh(...)");
|
|
|
|
out_node.resize(in_node.size());
|
|
for (int i = 0; i < out_node.size(); i++)
|
|
{
|
|
out_node[i].id = in_node[i].id;
|
|
out_node[i].x = in_node[i].x;
|
|
out_node[i].y = in_node[i].y;
|
|
|
|
out_node[i].time_ptr = in_time.get(i);
|
|
}
|
|
|
|
out_ele.resize(in_ele.size());
|
|
for (int i = 0; i < out_ele.size(); i++)
|
|
{
|
|
out_ele[i].init(in_ele[i].id,
|
|
out_node.get(in_ele[i].vert[0]->id),
|
|
out_node.get(in_ele[i].vert[1]->id),
|
|
out_node.get(in_ele[i].vert[2]->id),
|
|
in_slow.get(i));
|
|
}
|
|
return;
|
|
}
|
|
|
|
/**
|
|
* @brief 创建一个三维四面体网络的FMM正演网络
|
|
*
|
|
* @param[in] in_node 输入的三维顶点集
|
|
* @param[in] in_ele 输入的三维四面体集
|
|
* @param[in] in_time 输入的顶点顶点走时
|
|
* @param[in] in_slow 输入的四面体慢度集
|
|
* @param out_node 输出的三维FMM顶点集
|
|
* @param out_ele 输出的三维四面体集
|
|
*/
|
|
void gctl::create_fmm_mesh(const array<vertex3dc> &in_node, const array<tetrahedron> &in_ele, const array<double> &in_time,
|
|
const array<double> &in_slow, array<fmm_vertex3dc> &out_node, array<fmm_tetrahedron> &out_ele)
|
|
{
|
|
if (in_node.empty() || in_ele.empty() || in_time.empty() || in_slow.empty())
|
|
throw runtime_error("Invalid array sizes. From void gctl::create_fmm_mesh(...)");
|
|
|
|
if (in_node.size() != in_time.size())
|
|
throw runtime_error("Incompatible array sizes. From void gctl::create_fmm_mesh(...)");
|
|
|
|
out_node.resize(in_node.size());
|
|
for (int i = 0; i < out_node.size(); i++)
|
|
{
|
|
out_node[i].id = in_node[i].id;
|
|
out_node[i].x = in_node[i].x;
|
|
out_node[i].y = in_node[i].y;
|
|
out_node[i].z = in_node[i].z;
|
|
|
|
out_node[i].time_ptr = in_time.get(i);
|
|
}
|
|
|
|
out_ele.resize(in_ele.size());
|
|
for (int i = 0; i < out_ele.size(); i++)
|
|
{
|
|
out_ele[i].init(in_ele[i].id,
|
|
&out_node[in_ele[i].vert[0]->id],
|
|
&out_node[in_ele[i].vert[1]->id],
|
|
&out_node[in_ele[i].vert[2]->id],
|
|
&out_node[in_ele[i].vert[3]->id]);
|
|
out_ele[i].slow_ptr = in_slow.get(i);
|
|
}
|
|
return;
|
|
} |