zhangyiss/gctl_potential
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

tmp

Browse Source
This commit is contained in:
张壹 2025-01-09 18:57:08 +08:00
parent 2262b0269a
commit 6ba586004e
1 changed files with 93 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

99
lib/potential/gkernel_triangle2d.cpp
View File

@ -42,6 +42,8 @@ typedef void (*gkernel_triangle2d_ptr)(gctl::matrix<double> &out_kernel, const g
void gkernel_triangle2d_vz(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_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_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,
@ -65,6 +67,9 @@ void gctl::gkernel(matrix<double> &out_kernel, const array<triangle2d> &ele,
case Vz:
triangle_kernel = gkernel_triangle2d_vz;
break;
case Vx:
triangle_kernel = gkernel_triangle2d_vx;
break;
case Tzx:
triangle_kernel = gkernel_triangle2d_vzx;
break;
@ -92,6 +97,8 @@ typedef void (*gobser_triangle2d_ptr)(gctl::array<double> &out_obs, const gctl::
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);
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_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,
@ -116,6 +123,9 @@ void gctl::gobser(array<double> &out_obs, const array<triangle2d> &ele, const ar
case Vz:
triangle_obser = gobser_triangle2d_vz;
break;
case Vx:
triangle_obser = gobser_triangle2d_vx;
break;
case Tzx:
triangle_obser = gobser_triangle2d_vzx;
break;
@ -130,6 +140,7 @@ void gctl::gobser(array<double> &out_obs, const array<triangle2d> &ele, const ar
}
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_vzx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
double gkernel_triangle2d_vzz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
@ -174,7 +185,31 @@ void gkernel_triangle2d_vz(gctl::matrix<double> &out_kernel, const gctl::array<g
#pragma omp parallel for private (j) schedule(guided)
for (j = 0; j < e_size; j++)
{
out_kernel.at(i, j) = gkernel_triangle2d_vz_sig(ele.get(j), obsp.get(i));
out_kernel[i][j] = gkernel_triangle2d_vz_sig(ele.get(j), obsp.get(i));
}
}
return;
}
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)
{
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_vx");
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_vx_sig(ele.get(j), obsp.get(i));
}
}
return;
@ -198,7 +233,7 @@ void gkernel_triangle2d_vzx(gctl::matrix<double> &out_kernel, const gctl::array<
#pragma omp parallel for private (j) schedule(guided)
for (j = 0; j < e_size; j++)
{
out_kernel.at(i, j) = gkernel_triangle2d_vzx_sig(ele.get(j), obsp.get(i));
out_kernel[i][j] = gkernel_triangle2d_vzx_sig(ele.get(j), obsp.get(i));
}
}
return;
@ -222,7 +257,7 @@ void gkernel_triangle2d_vzz(gctl::matrix<double> &out_kernel, const gctl::array<
#pragma omp parallel for private (j) schedule(guided)
for (j = 0; j < e_size; j++)
{
out_kernel.at(i, j) = gkernel_triangle2d_vzz_sig(ele.get(j), obsp.get(i));
out_kernel[i][j] = gkernel_triangle2d_vzz_sig(ele.get(j), obsp.get(i));
}
}
return;
@ -246,7 +281,31 @@ void gobser_triangle2d_vz(gctl::array<double> &out_obs, const gctl::array<gctl::
#pragma omp parallel for private (i) schedule(guided)
for (i = 0; i < o_size; i++)
{
out_obs.at(i) += gkernel_triangle2d_vz_sig(ele.get(j), obsp.get(i)) * rho[j];
out_obs[i] += gkernel_triangle2d_vz_sig(ele.get(j), obsp.get(i)) * rho[j];
}
}
return;
}
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)
{
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_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) schedule(guided)
for (i = 0; i < o_size; i++)
{
out_obs[i] += gkernel_triangle2d_vx_sig(ele.get(j), obsp.get(i)) * rho[j];
}
}
return;
@ -270,7 +329,7 @@ void gobser_triangle2d_vzx(gctl::array<double> &out_obs, const gctl::array<gctl:
#pragma omp parallel for private (i) schedule(guided)
for (i = 0; i < o_size; i++)
{
out_obs.at(i) += gkernel_triangle2d_vzx_sig(ele.get(j), obsp.get(i)) * rho[j];
out_obs[i] += gkernel_triangle2d_vzx_sig(ele.get(j), obsp.get(i)) * rho[j];
}
}
return;
@ -294,7 +353,7 @@ void gobser_triangle2d_vzz(gctl::array<double> &out_obs, const gctl::array<gctl:
#pragma omp parallel for private (i) schedule(guided)
for (i = 0; i < o_size; i++)
{
out_obs.at(i) += gkernel_triangle2d_vzz_sig(ele.get(j), obsp.get(i)) * rho[j];
out_obs[i] += gkernel_triangle2d_vzz_sig(ele.get(j), obsp.get(i)) * rho[j];
}
}
return;
@ -328,6 +387,34 @@ double gkernel_triangle2d_vz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_p
return -2.0e+8*GCTL_G0*sum;
}
double gkernel_triangle2d_vx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr)
{
double sum;
double A, B, C;
double Aa, Ab, Ba, Bb, Ca, Cb;
sum = 0.0;
for (int n = 0; n < 3; n++)
{
Aa = (ele_ptr->vert[n]->y - op_ptr->y)*(ele_ptr->vert[(n+1)%3]->x - op_ptr->x) -
(ele_ptr->vert[(n+1)%3]->y - op_ptr->y)*(ele_ptr->vert[n]->x - op_ptr->x);
Ab = pow(ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y,2) +
pow(ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x,2);
A = Aa/Ab;
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);
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));
sum += A*(B+C);
}
return -2.0e+8*GCTL_G0*sum;
}
double gkernel_triangle2d_vzx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr)
{
double sum;