/********************************************************
* ██████╗ ██████╗████████╗██╗
* ██╔════╝ ██╔════╝╚══██╔══╝██║
* ██║ ███╗██║ ██║ ██║
* ██║ ██║██║ ██║ ██║
* ╚██████╔╝╚██████╗ ██║ ███████╗
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
* 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 "fmm_kernel.h"
#include "algorithm"
/**
* @brief 计算三角中第三个点的波前面走时
*
* 当三角形中的两个顶点波前面走时已知时可采用此公式计算第三个点的波前面走时。
* 假设三角形abc，其中顶点a与b的波前面走时已知（b的走时值大于a），求点c的波前面走时。
*
* @warning 此函数不应该暴露接口，不能在函数库中直接被调用。
*
* @param[in] area 三角形的面积
* @param[in] length_ab 边ab的长度
* @param[in] length_ac 边ac的长度
* @param[in] length_bc 边bc的长度
* @param[in] b_time 顶点b的波前面走时值
* @param[in] a_time 顶点a的波前面走时值
* @param[in] slow 三角形内的慢度值（波速的倒数）
* @param bk_grad 三角形内的走时相对于慢度的梯度
* @param. bk_ksi 三角形内的走时相对于慢度的梯度的系数
*
* @return 顶点c的波前面走时值
*/
double traveltime_triangle(double area, double length_ab, double length_ac,
double length_bc, double a_time, double b_time, double slow, double& bk_grad,
double& bk_ksi)
{
///< 沿边ac或bc传播的走时值
double direct_time;
double direct_ac_time = a_time + slow*length_ac;
double direct_bc_time = b_time + slow*length_bc;
if (direct_ac_time <= direct_bc_time)
{
direct_time = direct_ac_time; bk_grad = length_ac; bk_ksi = 0.0;
}
else
{
direct_time = direct_bc_time; bk_grad = length_bc; bk_ksi = 1.0;
}
double u = b_time - a_time;
double w = slow*slow*length_ab*length_ab - u*u;
if (w >= 0.0)
{
w = sqrt(w);
double rho0 = 2.0*area/length_ab;
double ksi0 = sqrt(length_ac*length_ac - rho0*rho0)/length_ab;
double ksi = ksi0 - u*rho0/(length_ab*w);
if (ksi > 0 && ksi < 1.0)
{
double trial_time = a_time + u*ksi0 + w*rho0/length_ab;
if (trial_time <= direct_time)
{
bk_grad = slow*length_ab*rho0/w;
bk_ksi = ksi;
return trial_time;
}
}
}
return direct_time;
}
/***************************************************/
/* */
/* 二维三角形网络的FMM算法 */
/* */
/***************************************************/
/**
* @brief 更新一个顶点的波前面走时值
*
* 循环顶点周围的所有三角形，若某个三角形中有两个顶点（非当前顶点）的走时已固定，
* 则可以利用已有走时信息计算当前点的走时值。如果新计算的走时值小于已有值，则替换
* 原有值。
*
* @param vert_ptr 顶点的指针
*/
void update_vertex_traveltime(gctl::fmm_vertex2dc* vert_ptr)
{
int j;
double time1, time2;
double back_time, back_grad, back_ksi;
// loop over all neighboring triangles. If any one have two fixed vertice, update time
for (int i = 0; i < vert_ptr->e_neigh.size(); i++)
{
// only proceed if the element is not used for calculation before.
if (vert_ptr->e_neigh[i]->tag == 0)
{
j = vert_ptr->ve_order[i];
if (vert_ptr->e_neigh[i]->vert[(j+1)%3]->tag == 2 &&
vert_ptr->e_neigh[i]->vert[(j+2)%3]->tag == 2)
{
time1 = *vert_ptr->e_neigh[i]->vert[(j+1)%3]->time_ptr;
time2 = *vert_ptr->e_neigh[i]->vert[(j+2)%3]->time_ptr;
if (time1 <= time2)
{
back_time = traveltime_triangle(vert_ptr->e_neigh[i]->area, vert_ptr->e_neigh[i]->edge_len[j],
vert_ptr->e_neigh[i]->edge_len[(j+2)%3], vert_ptr->e_neigh[i]->edge_len[(j+1)%3],
time1, time2, *vert_ptr->e_neigh[i]->slow_ptr, back_grad, back_ksi);
if (back_time < *vert_ptr->time_ptr)
{
*vert_ptr->time_ptr = back_time;
vert_ptr->jn_grad = back_grad;
vert_ptr->jn_ksi = back_ksi;
vert_ptr->jn_id1 = vert_ptr->e_neigh[i]->vert[(j+1)%3]->id;
vert_ptr->jn_id2 = vert_ptr->e_neigh[i]->vert[(j+2)%3]->id;
vert_ptr->jn_ele = vert_ptr->e_neigh[i]->id;
}
}
else
{
back_time = traveltime_triangle(vert_ptr->e_neigh[i]->area, vert_ptr->e_neigh[i]->edge_len[j],
vert_ptr->e_neigh[i]->edge_len[(j+1)%3], vert_ptr->e_neigh[i]->edge_len[(j+2)%3],
time2, time1, *vert_ptr->e_neigh[i]->slow_ptr, back_grad, back_ksi);
if (back_time < *vert_ptr->time_ptr)
{
*vert_ptr->time_ptr = back_time;
vert_ptr->jn_grad = back_grad;
vert_ptr->jn_ksi = back_ksi;
vert_ptr->jn_id1 = vert_ptr->e_neigh[i]->vert[(j+2)%3]->id;
vert_ptr->jn_id2 = vert_ptr->e_neigh[i]->vert[(j+1)%3]->id;
vert_ptr->jn_ele = vert_ptr->e_neigh[i]->id;
}
}
vert_ptr->e_neigh[i]->tag = 1;
}
}
}
return;
}
void gctl::fmm2d_forward_triangle(array *fmm_nodes,
array *fmm_triangles,
std::vector *front_node,
std::vector *record_node,
std::vector *end_node,
sparray2d *jn,
array *temp_jn)
{
// 检查运行条件至少有两个已知时间点且属于同一个三角形（某一个顶点的邻居也是确定点）
bool sec_break = false;
for (int i = 0; i < fmm_nodes->size(); i++)
{
if (fmm_nodes->at(i).tag == 2)
{
for (int j = 0; j < fmm_nodes->at(i).v_neigh.size(); j++)
{
// 条件满足
if (fmm_nodes->at(i).v_neigh[j]->tag == 2)
{
sec_break = true;
break;
}
}
}
if (sec_break) break;
}
// 条件未满足
if (!sec_break)
{
std::string err_str = "The FMM Fail to start. From void gctl::fmm2d_forward_triangle(...)";
throw runtime_error(err_str);
}
// 新建波前面
if (!front_node->empty()) front_node->clear();
for (int i = 0; i < fmm_nodes->size(); i++)
{
if (fmm_nodes->at(i).tag == 2)
{
for (int j = 0; j < fmm_nodes->at(i).v_neigh.size(); j++)
{
if (fmm_nodes->at(i).v_neigh[j]->tag == 0)
{
fmm_nodes->at(i).v_neigh[j]->tag = 1;
front_node->push_back(fmm_nodes->at(i).v_neigh[j]);
}
}
}
}
// 如果有三角形的三个顶点均已经固定 则将其设置为已计算状态
for (int i = 0; i < fmm_triangles->size(); i++)
{
if (fmm_triangles->at(i).vert[0]->tag == 2 &&
fmm_triangles->at(i).vert[1]->tag == 2 &&
fmm_triangles->at(i).vert[2]->tag == 2)
{
fmm_triangles->at(i).tag = 1;
}
}
//calculate trial time for all close nodes
for (int i = 0; i < front_node->size(); i++)
update_vertex_traveltime(front_node->at(i));
if (!record_node->empty()) record_node->clear();
//marching forward and updating the close nodes set
while (!front_node->empty())
{
// heap sort close nodes pointers to put the node first that has smallest time
//heap_sort_fmm2d(*front_node);
std::sort(front_node->begin(), front_node->end(), [](fmm_vertex2dc *a, fmm_vertex2dc *b)->bool{
if (*a->time_ptr < *b->time_ptr) return true; return false;
});
//change the first node's tag to 2 and update it's neighbor's time if it is not active
front_node->at(0)->tag = 2;
//记录节点的添加顺序
record_node->push_back(front_node->at(0));
for (int i = 0; i < front_node->at(0)->v_neigh.size(); i++)
{
if (front_node->at(0)->v_neigh[i]->tag == 0)
{
front_node->at(0)->v_neigh[i]->tag = 1;
update_vertex_traveltime(front_node->at(0)->v_neigh[i]);
front_node->push_back(front_node->at(0)->v_neigh[i]);
}
else if (front_node->at(0)->v_neigh[i]->tag == 1)
{
update_vertex_traveltime(front_node->at(0)->v_neigh[i]);
}
}
front_node->erase(front_node->begin());
// if all end_note are calculated. Stop the process
// otherwise, calculate all nodes' time.
if (end_node != nullptr)
{
std::sort(end_node->begin(), end_node->end(), [](fmm_vertex2dc *a, fmm_vertex2dc *b) -> bool{
if (a->tag < b->tag) return true; return false;
});
if (end_node->at(0)->tag == 2)
{
break;
}
}
}
// 节点走时对元素慢度的导数（雅各布矩阵，稀疏的）若无外部临时向量 则新建内部向量
if (jn != nullptr && temp_jn != nullptr)
{
//从节点记录的第0号位置开始计算JN
for (int i = 0; i < record_node->size(); i++)
{
for (int j = 0; j < fmm_triangles->size(); j++)
temp_jn->at(j) = 0.0;
if (record_node->at(i)->jn_ksi != 1.0)
{
jn->at(record_node->at(i)->jn_id1)->export_dense(*temp_jn, 1.0-record_node->at(i)->jn_ksi, AppendVal);
}
if (record_node->at(i)->jn_ksi != 0.0)
{
jn->at(record_node->at(i)->jn_id2)->export_dense(*temp_jn, record_node->at(i)->jn_ksi, AppendVal);
}
temp_jn->at(record_node->at(i)->jn_ele) += record_node->at(i)->jn_grad;
jn->at(record_node->at(i)->id)->malloc(*temp_jn, 0.0, GCTL_ZERO);
}
}
return;
}
void gctl::source2receiver_direct(array *fmm_nodes, array *fmm_triangles,
seis_point2d_tri *in_source, seis_point2d_tri *in_receiver, array *rece_ele_grad,
std::vector *close_node_ptr, std::vector *record_node_ptr,
sparray2d *jn_ptr, array *tmp_jn_ptr)
{
// 初始化波前面记录向量 若无外部临时向量 则新建内部向量
std::vector *front_node_ptr = nullptr;
if (close_node_ptr != nullptr)
{
front_node_ptr = close_node_ptr;
if (!front_node_ptr->empty()) front_node_ptr->clear();
}
else front_node_ptr = new std::vector;
// reserve space
front_node_ptr->reserve(fmm_nodes->size());
// 初始化节点旅行记录向量 若无外部临时向量 则新建内部向量
std::vector *travel_record_ptr = nullptr;
if (record_node_ptr != nullptr)
{
travel_record_ptr = record_node_ptr;
if (!travel_record_ptr->empty()) travel_record_ptr->clear();
}
else travel_record_ptr = new std::vector;
// reserve space
travel_record_ptr->reserve(fmm_nodes->size());
sparray2d *JN = nullptr;
// 节点走时对元素慢度的导数（雅各布矩阵，稀疏的）若无外部临时向量 则新建内部向量
if (rece_ele_grad != nullptr)
{
if (jn_ptr != nullptr)
{
JN = jn_ptr;
for (int i = 0; i < JN->size(); i++)
{
if (!JN->at(i)->empty()) JN->at(i)->clear();
}
}
else JN = new sparray2d(fmm_nodes->size(), fmm_triangles->size(), 0.0);
}
array *temp_JN = nullptr;
// 若无外部临时向量 则新建内部向量
if (rece_ele_grad != nullptr)
{
if (tmp_jn_ptr != nullptr)
temp_JN = tmp_jn_ptr;
else temp_JN = new array(fmm_triangles->size());
}
// 重置所有标记与时间 从激发点开始计算
for (int i = 0; i < fmm_nodes->size(); i++)
{
fmm_nodes->at(i).tag = 0; *fmm_nodes->at(i).time_ptr = GCTL_BDL_MAX;
}
for (int i = 0; i < fmm_triangles->size(); i++)
{
fmm_triangles->at(i).tag = 0;
}
// find host elements for source and receiver
in_source->find_host_element(fmm_triangles->get(), fmm_triangles->size());
in_receiver->find_host_element(fmm_triangles->get(), fmm_triangles->size());
// assign initial tags for elements
in_source->host_ele->tag = 1;
for (int i = 0; i < 3; i++)
{
in_source->host_ele->vert[i]->tag = 2;
*in_source->host_ele->vert[i]->time_ptr = *in_source->host_ele->slow_ptr *
distance(*in_source->host_ele->vert[i], *in_source);
}
double temp_val;
if (rece_ele_grad != nullptr)
{
for (int i = 0; i < 3; i++)
{
//初始化前三个梯度值
//对元素的慢度求梯度即为源到顶点的距离
temp_val= distance(*in_source->host_ele->vert[i], *in_source);
JN->at(in_source->host_ele->vert[i]->id)->set(in_source->host_ele->id, temp_val);
}
}
std::vector rece_nodes;
for (int i = 0; i < 3; i++)
{
rece_nodes.push_back(in_receiver->host_ele->vert[i]);
}
// Calculate using FMM
fmm2d_forward_triangle(fmm_nodes, fmm_triangles, front_node_ptr, travel_record_ptr,
&rece_nodes, JN, temp_JN);
double r[3], w_sum;
for (int i = 0; i < 3; i++)
{
r[i] = distance(*in_receiver->host_ele->vert[i],
*in_receiver) + GCTL_ZERO;
}
w_sum = 1.0/r[0] + 1.0/r[1] + 1.0/r[2];
in_receiver->time = *in_receiver->host_ele->vert[0]->time_ptr/(r[0]*w_sum) +
*in_receiver->host_ele->vert[1]->time_ptr/(r[1]*w_sum) +
*in_receiver->host_ele->vert[2]->time_ptr/(r[2]*w_sum);
if (rece_ele_grad != nullptr)
{
//计算接收点的JN
for (int i = 0; i < fmm_triangles->size(); i++)
rece_ele_grad->at(i) = 0.0;
for (int i = 0; i < 3; i++)
{
JN->at(in_receiver->host_ele->vert[i]->id)->export_dense(*rece_ele_grad, 1.0/(r[i]*w_sum), AppendVal);
}
}
if (record_node_ptr == nullptr) delete travel_record_ptr;
if (close_node_ptr == nullptr) delete front_node_ptr;
if (jn_ptr == nullptr && JN != nullptr) delete JN;
if (tmp_jn_ptr == nullptr && temp_JN != nullptr) delete temp_JN;
return;
}
void gctl::source2receiver_reflect(array *fmm_nodes, array *fmm_triangles,
seis_point2d_tri *in_source, seis_point2d_tri *in_receiver, std::vector *reflect_node,
array *rece_ele_grad, std::vector *close_node_ptr,
std::vector *record_node_ptr, sparray2d *jn_ptr,
array *tmp_jn_ptr, array *reflect_index)
{
// 初始化波前面记录向量 若无外部临时向量 则新建内部向量
std::vector *front_node_ptr = nullptr;
if (close_node_ptr != nullptr)
{
front_node_ptr = close_node_ptr;
if (!front_node_ptr->empty()) front_node_ptr->clear();
}
else front_node_ptr = new std::vector;
// reserve space
front_node_ptr->reserve(fmm_nodes->size());
// 初始化节点旅行记录向量 若无外部临时向量 则新建内部向量
std::vector *travel_record_ptr = nullptr;
if (record_node_ptr != nullptr)
{
travel_record_ptr = record_node_ptr;
if (!travel_record_ptr->empty()) travel_record_ptr->clear();
}
else travel_record_ptr = new std::vector;
// reserve space
travel_record_ptr->reserve(fmm_nodes->size());
sparray2d *JN = nullptr;
// 节点走时对元素慢度的导数（雅各布矩阵，稀疏的）若无外部临时向量 则新建内部向量
if (rece_ele_grad != nullptr)
{
if (jn_ptr != nullptr)
{
JN = jn_ptr;
for (int i = 0; i < JN->size(); i++)
{
if (!JN->at(i)->empty()) JN->at(i)->clear();
}
}
else JN = new sparray2d(fmm_nodes->size(), fmm_triangles->size(), 0.0);
}
array *temp_JN = nullptr;
// 若无外部临时向量 则新建内部向量
if (rece_ele_grad != nullptr)
{
if (tmp_jn_ptr != nullptr)
temp_JN = tmp_jn_ptr;
else temp_JN = new array(fmm_triangles->size());
}
array *non_reflect_index = nullptr;
if (reflect_index != nullptr)
non_reflect_index = reflect_index;
else non_reflect_index = new array(fmm_nodes->size());
// 重置所有标记与时间 从激发点开始计算
for (int i = 0; i < fmm_nodes->size(); i++)
{
fmm_nodes->at(i).tag = 0; *fmm_nodes->at(i).time_ptr = GCTL_BDL_MAX;
non_reflect_index->at(i) = true;
}
for (int i = 0; i < fmm_triangles->size(); i++)
{
fmm_triangles->at(i).tag = 0;
}
// find host elements for source and receiver
in_source->find_host_element(fmm_triangles->get(), fmm_triangles->size());
in_receiver->find_host_element(fmm_triangles->get(), fmm_triangles->size());
// assign initial tags for elements
in_source->host_ele->tag = 1;
for (int i = 0; i < 3; i++)
{
in_source->host_ele->vert[i]->tag = 2;
*in_source->host_ele->vert[i]->time_ptr = *in_source->host_ele->slow_ptr *
distance(*in_source->host_ele->vert[i], *in_source);
}
double temp_val;
if (rece_ele_grad != nullptr)
{
for (int i = 0; i < 3; i++)
{
//初始化前三个梯度值
//对元素的慢度求梯度即为源到顶点的距离
temp_val= distance(*in_source->host_ele->vert[i], *in_source);
JN->at(in_source->host_ele->vert[i]->id)->set(in_source->host_ele->id, temp_val);
}
}
// 第一次计算的接收点为反射面上的点 记录那些点是非反射界面上的点
std::vector rece_nodes;
for (int i = 0; i < reflect_node->size(); i++)
{
non_reflect_index->at(reflect_node->at(i)->id) = false;
rece_nodes.push_back(reflect_node->at(i));
}
// Calculate using FMM
fmm2d_forward_triangle(fmm_nodes, fmm_triangles, front_node_ptr, travel_record_ptr,
&rece_nodes, JN, temp_JN);
// 重置出反射点之外的所有点
for (int i = 0; i < fmm_nodes->size(); i++)
{
if (non_reflect_index->at(i))
{
fmm_nodes->at(i).tag = 0; *fmm_nodes->at(i).time_ptr = GCTL_BDL_MAX;
JN->at(i)->clear();
}
}
for (int i = 0; i < fmm_triangles->size(); i++)
{
fmm_triangles->at(i).tag = 0;
}
// 第二次计算为接收点
for (int i = 0; i < 3; i++)
{
rece_nodes.push_back(in_receiver->host_ele->vert[i]);
}
if (!front_node_ptr->empty()) front_node_ptr->clear();
if (!travel_record_ptr->empty()) travel_record_ptr->clear();
// Calculate using FMM
fmm2d_forward_triangle(fmm_nodes, fmm_triangles, front_node_ptr, travel_record_ptr,
&rece_nodes, JN, temp_JN);
double r[3], w_sum;
for (int i = 0; i < 3; i++)
{
r[i] = distance(*in_receiver->host_ele->vert[i],
*in_receiver) + GCTL_ZERO;
}
w_sum = 1.0/r[0] + 1.0/r[1] + 1.0/r[2];
in_receiver->time = *in_receiver->host_ele->vert[0]->time_ptr/(r[0]*w_sum) +
*in_receiver->host_ele->vert[1]->time_ptr/(r[1]*w_sum) +
*in_receiver->host_ele->vert[2]->time_ptr/(r[2]*w_sum);
if (rece_ele_grad != nullptr)
{
//计算接收点的JN
for (int i = 0; i < fmm_triangles->size(); i++)
rece_ele_grad->at(i) = 0.0;
for (int i = 0; i < 3; i++)
{
JN->at(in_receiver->host_ele->vert[i]->id)->export_dense(*rece_ele_grad, 1.0/(r[i]*w_sum), AppendVal);
}
}
if (record_node_ptr == nullptr) delete travel_record_ptr;
if (close_node_ptr == nullptr) delete front_node_ptr;
if (jn_ptr == nullptr && JN != nullptr) delete JN;
if (tmp_jn_ptr == nullptr && temp_JN != nullptr) delete temp_JN;
if (reflect_index == nullptr) delete non_reflect_index;
return;
}
/***************************************************/
/* */
/* 三维四面体网络的FMM算法 */
/* */
/***************************************************/
double traveltime_tetrahedron(gctl::fmm_tetrahedron* ele_ptr, int d_order, int a_order,
int b_order, int c_order, double a_time, double b_time, double c_time,
double &bk_grad, double &bk_ksi, double &bk_zeta)
{
double length_ab, length_ac, length_bc, length_ad, length_bd, length_cd;
double area_abc, area_adb, area_bdc, area_cda;
length_ab = ele_ptr->edge_len[ele_ptr->vert2edge_id(a_order, b_order)];
length_ac = ele_ptr->edge_len[ele_ptr->vert2edge_id(a_order, c_order)];
length_bc = ele_ptr->edge_len[ele_ptr->vert2edge_id(b_order, c_order)];
length_ad = ele_ptr->edge_len[ele_ptr->vert2edge_id(a_order, d_order)];
length_bd = ele_ptr->edge_len[ele_ptr->vert2edge_id(b_order, d_order)];
length_cd = ele_ptr->edge_len[ele_ptr->vert2edge_id(c_order, d_order)];
area_abc = ele_ptr->face_area[ele_ptr->vert2face_id(a_order, b_order, c_order)];
area_adb = ele_ptr->face_area[ele_ptr->vert2face_id(a_order, d_order, b_order)];
area_bdc = ele_ptr->face_area[ele_ptr->vert2face_id(b_order, d_order, c_order)];
area_cda = ele_ptr->face_area[ele_ptr->vert2face_id(c_order, d_order, a_order)];
bk_zeta = 0.0;
double slow = *ele_ptr->slow_ptr;
double direct_time, tmp_time, tmp_bkg, tmp_bkk;
direct_time = traveltime_triangle(area_adb, length_ab, length_ad, length_bd,
a_time, b_time, slow, bk_grad, bk_ksi);
tmp_time = traveltime_triangle(area_cda, length_ac, length_ad, length_cd,
a_time, c_time, slow, tmp_bkg, tmp_bkk);
if (tmp_time < direct_time)
{
direct_time = tmp_time;
bk_grad = tmp_bkg;
bk_ksi = tmp_bkk;
}
tmp_time = traveltime_triangle(area_bdc, length_bc, length_bd, length_cd,
b_time, c_time, slow, tmp_bkg, tmp_bkk);
if (tmp_time < direct_time)
{
direct_time = tmp_time;
bk_grad = tmp_bkg;
bk_ksi = tmp_bkk;
}
double phi = 2.0*area_abc;
double u = b_time - a_time;
double v = c_time - a_time;
double dd = 0.5*(length_ab*length_ab + length_ac*length_ac - length_bc*length_bc);
double w = slow*slow*phi*phi - u*u*length_ac*length_ac - v*v*length_ab*length_ab + 2.0*u*v*dd;
if (w >= 0.0)
{
w = sqrt(w);
double rho0 = 3.0*ele_ptr->volume/area_abc;
double ksi0 = sqrt(4.0*area_cda*area_cda/(length_ac*length_ac) - rho0*rho0)/(2.0*area_abc/length_ac);
double zeta0 = sqrt(4.0*area_adb*area_adb/(length_ab*length_ab) - rho0*rho0)/(2.0*area_abc/length_ab);
double ksi = ksi0 - fabs(u*length_ac*length_ac - v*dd)*rho0/(phi*w);
double zeta = zeta0 - fabs(v*length_ab*length_ab - u*dd)*rho0/(phi*w);
if (ksi > 0 && ksi < 1 && zeta > 0 && zeta < 1 && ksi+zeta > 0 && ksi+zeta < 1)
{
double trial_time = a_time + u*ksi0 + v*zeta0 + w*rho0/phi;
if (trial_time > c_time && trial_time < direct_time) // c_time是最大的 a b c 里面最大的
{
bk_grad = slow*phi*rho0/w;
bk_ksi = ksi;
bk_zeta = zeta;
return trial_time;
}
}
}
return direct_time;
}
/**
* @brief 更新一个顶点的波前面走时值
*
* 循环顶点周围的所有四面体，若某个四面体中有三个顶点（非当前顶点）的走时已固定，
* 则可以利用已有走时信息计算当前点的走时值。如果新计算的走时值小于已有值，则替换
* 原有值。
*
* @param vert_ptr 顶点的指针
*/
void update_vertex_traveltime(gctl::fmm_vertex3dc* vert_ptr)
{
int j;
double time1, time2, time3;
double back_time, back_grad, back_ksi, back_zeta;
// loop over all neighboring triangles. If any one have two fixed vertice, update time
for (int i = 0; i < vert_ptr->e_neigh.size(); i++)
{
// only proceed if the element is not used for calculation before.
if (vert_ptr->e_neigh[i]->tag == 0)
{
j = vert_ptr->ve_order[i];
if (vert_ptr->e_neigh[i]->vert[(j+1)%4]->tag == 2 &&
vert_ptr->e_neigh[i]->vert[(j+2)%4]->tag == 2 &&
vert_ptr->e_neigh[i]->vert[(j+3)%4]->tag == 2)
{
time1 = *vert_ptr->e_neigh[i]->vert[(j+1)%4]->time_ptr;
time2 = *vert_ptr->e_neigh[i]->vert[(j+2)%4]->time_ptr;
time3 = *vert_ptr->e_neigh[i]->vert[(j+3)%4]->time_ptr;
if (time1 <= time2 && time2 <= time3) //123
{
back_time = traveltime_tetrahedron(vert_ptr->e_neigh[i], j, (j+1)%4,
(j+2)%4, (j+3)%4, time1, time2, time3, back_grad, back_ksi, back_zeta);
if (back_time < *vert_ptr->time_ptr)
{
*vert_ptr->time_ptr = back_time;
vert_ptr->jn_grad = back_grad;
vert_ptr->jn_ksi = back_ksi;
vert_ptr->jn_zeta = back_zeta;
vert_ptr->jn_id1 = vert_ptr->e_neigh[i]->vert[(j+1)%4]->id;
vert_ptr->jn_id2 = vert_ptr->e_neigh[i]->vert[(j+2)%4]->id;
vert_ptr->jn_id3 = vert_ptr->e_neigh[i]->vert[(j+3)%4]->id;
vert_ptr->jn_ele = vert_ptr->e_neigh[i]->id;
}
}
else if (time1 <= time3 && time3 <= time2) //132
{
back_time = traveltime_tetrahedron(vert_ptr->e_neigh[i], j, (j+1)%4,
(j+3)%4, (j+2)%4, time1, time3, time2, back_grad, back_ksi, back_zeta);
if (back_time < *vert_ptr->time_ptr)
{
*vert_ptr->time_ptr = back_time;
vert_ptr->jn_grad = back_grad;
vert_ptr->jn_ksi = back_ksi;
vert_ptr->jn_zeta = back_zeta;
vert_ptr->jn_id1 = vert_ptr->e_neigh[i]->vert[(j+1)%4]->id;
vert_ptr->jn_id2 = vert_ptr->e_neigh[i]->vert[(j+3)%4]->id;
vert_ptr->jn_id3 = vert_ptr->e_neigh[i]->vert[(j+2)%4]->id;
vert_ptr->jn_ele = vert_ptr->e_neigh[i]->id;
}
}
else if (time2 <= time1 && time1 <= time3) //213
{
back_time = traveltime_tetrahedron(vert_ptr->e_neigh[i], j, (j+2)%4,
(j+1)%4, (j+3)%4, time2, time1, time3, back_grad, back_ksi, back_zeta);
if (back_time < *vert_ptr->time_ptr)
{
*vert_ptr->time_ptr = back_time;
vert_ptr->jn_grad = back_grad;
vert_ptr->jn_ksi = back_ksi;
vert_ptr->jn_zeta = back_zeta;
vert_ptr->jn_id1 = vert_ptr->e_neigh[i]->vert[(j+2)%4]->id;
vert_ptr->jn_id2 = vert_ptr->e_neigh[i]->vert[(j+1)%4]->id;
vert_ptr->jn_id3 = vert_ptr->e_neigh[i]->vert[(j+3)%4]->id;
vert_ptr->jn_ele = vert_ptr->e_neigh[i]->id;
}
}
else if (time2 <= time3 && time3 <= time1) //231
{
back_time = traveltime_tetrahedron(vert_ptr->e_neigh[i], j, (j+2)%4,
(j+3)%4, (j+1)%4, time2, time3, time1, back_grad, back_ksi, back_zeta);
if (back_time < *vert_ptr->time_ptr)
{
*vert_ptr->time_ptr = back_time;
vert_ptr->jn_grad = back_grad;
vert_ptr->jn_ksi = back_ksi;
vert_ptr->jn_zeta = back_zeta;
vert_ptr->jn_id1 = vert_ptr->e_neigh[i]->vert[(j+2)%4]->id;
vert_ptr->jn_id2 = vert_ptr->e_neigh[i]->vert[(j+3)%4]->id;
vert_ptr->jn_id3 = vert_ptr->e_neigh[i]->vert[(j+1)%4]->id;
vert_ptr->jn_ele = vert_ptr->e_neigh[i]->id;
}
}
else if (time3 <= time1 && time1 <= time2) //312
{
back_time = traveltime_tetrahedron(vert_ptr->e_neigh[i], j, (j+3)%4,
(j+1)%4, (j+2)%4, time3, time1, time2, back_grad, back_ksi, back_zeta);
if (back_time < *vert_ptr->time_ptr)
{
*vert_ptr->time_ptr = back_time;
vert_ptr->jn_grad = back_grad;
vert_ptr->jn_ksi = back_ksi;
vert_ptr->jn_zeta = back_zeta;
vert_ptr->jn_id1 = vert_ptr->e_neigh[i]->vert[(j+3)%4]->id;
vert_ptr->jn_id2 = vert_ptr->e_neigh[i]->vert[(j+1)%4]->id;
vert_ptr->jn_id3 = vert_ptr->e_neigh[i]->vert[(j+2)%4]->id;
vert_ptr->jn_ele = vert_ptr->e_neigh[i]->id;
}
}
else //321
{
back_time = traveltime_tetrahedron(vert_ptr->e_neigh[i], j, (j+3)%4,
(j+2)%4, (j+1)%4, time3, time2, time1, back_grad, back_ksi, back_zeta);
if (back_time < *vert_ptr->time_ptr)
{
*vert_ptr->time_ptr = back_time;
vert_ptr->jn_grad = back_grad;
vert_ptr->jn_ksi = back_ksi;
vert_ptr->jn_zeta = back_zeta;
vert_ptr->jn_id1 = vert_ptr->e_neigh[i]->vert[(j+3)%4]->id;
vert_ptr->jn_id2 = vert_ptr->e_neigh[i]->vert[(j+2)%4]->id;
vert_ptr->jn_id3 = vert_ptr->e_neigh[i]->vert[(j+1)%4]->id;
vert_ptr->jn_ele = vert_ptr->e_neigh[i]->id;
}
}
vert_ptr->e_neigh[i]->tag = 1; // set tag to 1 after calculation
}
}
}
return;
}
void gctl::fmm3d_forward_tetrahedron(array *fmm_nodes,
array *fmm_tetrahedrons,
std::vector *front_node,
std::vector *record_node,
std::vector *end_node,
sparray2d *jn,
array *temp_jn)
{
// 检查运行条件 至少有三个已知时间点且属于同一个四面体
bool run_ok = false;
int valid_count;
for (int i = 0; i < fmm_tetrahedrons->size(); i++)
{
valid_count = 0;
for (int j = 0; j < 4; j++)
{
if (fmm_tetrahedrons->at(i).vert[j]->tag == 2)
{
valid_count++;
}
}
if (valid_count >= 3)
{
run_ok = true;
break;
}
}
// 条件未满足
if (!run_ok)
{
throw runtime_error("FMM fail to start. From fmm3d_forward_tetrahedron(...)");
}
if (!front_node->empty()) front_node->clear();
// 新建波前面
for (int i = 0; i < fmm_nodes->size(); i++)
{
if (fmm_nodes->at(i).tag == 2)
{
for (int j = 0; j < fmm_nodes->at(i).v_neigh.size(); j++)
{
if (fmm_nodes->at(i).v_neigh[j]->tag == 0)
{
fmm_nodes->at(i).v_neigh[j]->tag = 1;
front_node->push_back(fmm_nodes->at(i).v_neigh[j]);
}
}
}
}
// 如果有四面体的四个顶点均已经固定 则将其设置为已计算状态
for (int i = 0; i < fmm_tetrahedrons->size(); i++)
{
if (fmm_tetrahedrons->at(i).vert[0]->tag == 2 &&
fmm_tetrahedrons->at(i).vert[1]->tag == 2 &&
fmm_tetrahedrons->at(i).vert[2]->tag == 2 &&
fmm_tetrahedrons->at(i).vert[3]->tag == 2)
{
fmm_tetrahedrons->at(i).tag = 1;
}
}
//calculate trial time for all close nodes
for (int i = 0; i < front_node->size(); i++)
update_vertex_traveltime(front_node->at(i));
if (!record_node->empty()) record_node->clear();
//marching forward and updating the close nodes set
while (!front_node->empty())
{
// heap sort close nodes pointers to put the node first that has smallest time
std::sort(front_node->begin(), front_node->end(), [](fmm_vertex3dc *a, fmm_vertex3dc *b)->bool{
if (*a->time_ptr < *b->time_ptr) return true; return false;
});
//change the first node's tag to 2 and update it's neighbor's time if it is not active
front_node->at(0)->tag = 2;
//记录节点的添加顺序
record_node->push_back(front_node->at(0));
for (int i = 0; i < front_node->at(0)->v_neigh.size(); i++)
{
if (front_node->at(0)->v_neigh[i]->tag == 0)
{
front_node->at(0)->v_neigh[i]->tag = 1;
update_vertex_traveltime(front_node->at(0)->v_neigh[i]);
front_node->push_back(front_node->at(0)->v_neigh[i]);
}
else if (front_node->at(0)->v_neigh[i]->tag == 1)
{
update_vertex_traveltime(front_node->at(0)->v_neigh[i]);
}
}
front_node->erase(front_node->begin());
// if all end_note are calculated. Stop the process. Otherwise, calculate all nodes' time.
if (end_node != nullptr)
{
std::sort(end_node->begin(), end_node->end(), [](fmm_vertex3dc *a, fmm_vertex3dc *b)->bool{
if (a->tag < b->tag) return true; return false;
});
if (end_node->at(0)->tag == 2)
{
break;
}
}
}
// 节点走时对元素慢度的导数（雅各布矩阵，稀疏的）若无外部临时向量 则新建内部向量
if (jn != nullptr && temp_jn != nullptr)
{
//从节点记录的第0号位置开始计算JN
for (int i = 0; i < record_node->size(); i++)
{
for (int j = 0; j < fmm_tetrahedrons->size(); j++)
temp_jn->at(j) = 0.0;
if (record_node->at(i)->jn_ksi != 1.0)
{
jn->at(record_node->at(i)->jn_id1)->export_dense(*temp_jn, 1.0 - record_node->at(i)->jn_ksi - record_node->at(i)->jn_zeta, AppendVal);
}
if (record_node->at(i)->jn_ksi != 0.0)
{
jn->at(record_node->at(i)->jn_id2)->export_dense(*temp_jn, record_node->at(i)->jn_ksi, AppendVal);
}
if (record_node->at(i)->jn_zeta != 0.0)
{
jn->at(record_node->at(i)->jn_id3)->export_dense(*temp_jn, record_node->at(i)->jn_zeta, AppendVal);
}
temp_jn->at(record_node->at(i)->jn_ele) += record_node->at(i)->jn_grad;
jn->at(record_node->at(i)->id)->malloc(*temp_jn, 0.0, GCTL_ZERO);
}
}
return;
}
void gctl::source2receiver_direct(array *fmm_nodes, array *fmm_tetrahedrons,
seis_point3d_tet *in_source, seis_point3d_tet *in_receiver, array *rece_ele_grad,
std::vector *close_node_ptr, std::vector *record_node_ptr,
sparray2d *jn_ptr, array *tmp_jn_ptr)
{
// 初始化波前面记录向量 若无外部临时向量 则新建内部向量
std::vector *front_node_ptr = nullptr;
if (close_node_ptr != nullptr)
{
front_node_ptr = close_node_ptr;
if (!front_node_ptr->empty()) front_node_ptr->clear();
}
else front_node_ptr = new std::vector;
// reserve space
front_node_ptr->reserve(fmm_nodes->size());
// 初始化节点旅行记录向量 若无外部临时向量 则新建内部向量
std::vector *travel_record_ptr = nullptr;
if (record_node_ptr != nullptr)
{
travel_record_ptr = record_node_ptr;
if (!travel_record_ptr->empty()) travel_record_ptr->clear();
}
else travel_record_ptr = new std::vector;
// reserve space
travel_record_ptr->reserve(fmm_nodes->size());
sparray2d *JN = nullptr;
// 节点走时对元素慢度的导数（雅各布矩阵，稀疏的）若无外部临时向量 则新建内部向量
if (rece_ele_grad != nullptr)
{
if (jn_ptr != nullptr)
{
JN = jn_ptr;
for (int i = 0; i < JN->size(); i++)
{
if (!JN->at(i)->empty()) JN->at(i)->clear();
}
}
else JN = new sparray2d(fmm_nodes->size(), fmm_tetrahedrons->size(), 0.0);
}
array *temp_JN = nullptr;
// 若无外部临时向量 则新建内部向量
if (rece_ele_grad != nullptr)
{
if (tmp_jn_ptr != nullptr)
temp_JN = tmp_jn_ptr;
else temp_JN = new array(fmm_tetrahedrons->size());
}
// 重置所有标记与时间 从激发点开始计算
for (int i = 0; i < fmm_nodes->size(); i++)
{
fmm_nodes->at(i).tag = 0; *fmm_nodes->at(i).time_ptr = GCTL_BDL_MAX;
}
for (int i = 0; i < fmm_tetrahedrons->size(); i++)
{
fmm_tetrahedrons->at(i).tag = 0;
}
// find host elements for source and receiver
in_source->find_host_element(fmm_tetrahedrons->get(), fmm_tetrahedrons->size());
in_receiver->find_host_element(fmm_tetrahedrons->get(), fmm_tetrahedrons->size());
// assign initial tags for elements
in_source->host_ele->tag = 1;
for (int i = 0; i < 4; i++)
{
in_source->host_ele->vert[i]->tag = 2;
*in_source->host_ele->vert[i]->time_ptr = *in_source->host_ele->slow_ptr *
distance(*in_source->host_ele->vert[i], *in_source);
}
double temp_val;
if (rece_ele_grad != nullptr)
{
for (int i = 0; i < 4; i++)
{
//初始化前三个梯度值
//对元素的慢度求梯度即为源到顶点的距离
temp_val= distance(*in_source->host_ele->vert[i], *in_source);
JN->at(in_source->host_ele->vert[i]->id)->set(in_source->host_ele->id, temp_val);
}
}
std::vector rece_nodes;
for (int i = 0; i < 4; i++)
{
rece_nodes.push_back(in_receiver->host_ele->vert[i]);
}
// Calculate using FMM
fmm3d_forward_tetrahedron(fmm_nodes, fmm_tetrahedrons, front_node_ptr, travel_record_ptr,
&rece_nodes, JN, temp_JN);
double r[4], w_sum;
for (int i = 0; i < 4; i++)
{
r[i] = distance(*in_receiver->host_ele->vert[i],
*in_receiver) + GCTL_ZERO;
}
w_sum = 1.0/r[0] + 1.0/r[1] + 1.0/r[2] + 1.0/r[3];
in_receiver->time = *in_receiver->host_ele->vert[0]->time_ptr/(r[0]*w_sum) +
*in_receiver->host_ele->vert[1]->time_ptr/(r[1]*w_sum) +
*in_receiver->host_ele->vert[2]->time_ptr/(r[2]*w_sum) +
*in_receiver->host_ele->vert[3]->time_ptr/(r[3]*w_sum);
if (rece_ele_grad != nullptr)
{
//计算接收点的JN
for (int i = 0; i < fmm_tetrahedrons->size(); i++)
rece_ele_grad->at(i) = 0.0;
for (int i = 0; i < 4; i++)
{
JN->at(in_receiver->host_ele->vert[i]->id)->export_dense(*rece_ele_grad, 1.0/(r[i]*w_sum), AppendVal);
}
}
if (record_node_ptr == nullptr) delete travel_record_ptr;
if (close_node_ptr == nullptr) delete front_node_ptr;
if (jn_ptr == nullptr && JN != nullptr) delete JN;
if (tmp_jn_ptr == nullptr && temp_JN != nullptr) delete temp_JN;
return;
}