zhangyiss/gctl_potential
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
0988ff4a7f
gctl_potential/lib/potential/gkernel_block.cpp
Yi Zhang 042d310021 initial upload
2024-09-10 19:56:41 +08:00

530 lines
21 KiB
C++
Raw Blame History

/********************************************************
* ██████╗ ██████╗████████╗██╗
* ██╔════╝ ██╔════╝╚══██╔══╝██║
* ██║ ███╗██║ ██║ ██║
* ██║ ██║██║ ██║ ██║
* ╚██████╔╝╚██████╗ ██║ ███████╗
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
* 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 "gkernel_block.h"
#include "cmath"
typedef void (*gkernel_block_ptr)(gctl::matrix<double> &out_kernel, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, gctl::verbose_type_e verbose);
void gkernel_block_vz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, gctl::verbose_type_e verbose);
void gkernel_block_vzx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, gctl::verbose_type_e verbose);
void gkernel_block_vzy(gctl::matrix<double> &out_kernel, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, gctl::verbose_type_e verbose);
void gkernel_block_vzz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, gctl::verbose_type_e verbose);
void gctl::gkernel(matrix<double> &out_kernel, const array<block> &ele, const array<point3dc> &obsp,
gravitational_field_type_e comp_id, verbose_type_e verbose)
{
gkernel_block_ptr block_kernel;
switch (comp_id)
{
case Vz:
block_kernel = gkernel_block_vz;
break;
case Tzx:
block_kernel = gkernel_block_vzx;
break;
case Tzy:
block_kernel = gkernel_block_vzy;
break;
case Tzz:
block_kernel = gkernel_block_vzz;
break;
default:
block_kernel = gkernel_block_vz;
break;
}
block_kernel(out_kernel, ele, obsp, verbose);
return;
}
typedef void (*gobser_block_ptr)(gctl::array<double> &out_obs, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
void gobser_block_vz(gctl::array<double> &out_obs, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
void gobser_block_vzx(gctl::array<double> &out_obs, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
void gobser_block_vzy(gctl::array<double> &out_obs, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
void gobser_block_vzz(gctl::array<double> &out_obs, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verboses);
void gctl::gobser(array<double> &out_obs, const array<block> &ele, const array<point3dc> &obsp,
const array<double> &rho, gravitational_field_type_e comp_id, verbose_type_e verbose)
{
gobser_block_ptr block_obser;
switch (comp_id)
{
case Vz:
block_obser = gobser_block_vz;
break;
case Tzx:
block_obser = gobser_block_vzx;
break;
case Tzy:
block_obser = gobser_block_vzy;
break;
case Tzz:
block_obser = gobser_block_vzz;
break;
default:
block_obser = gobser_block_vz;
break;
}
block_obser(out_obs, ele, obsp, rho, verbose);
return;
}
// declare algorithm for individual element and observation point
double gkernel_block_vz_sig(gctl::block *ele_ptr, gctl::point3dc *op_ptr);
double gkernel_block_vzx_sig(gctl::block *ele_ptr, gctl::point3dc *op_ptr);
double gkernel_block_vzy_sig(gctl::block *ele_ptr, gctl::point3dc *op_ptr);
double gkernel_block_vzz_sig(gctl::block *ele_ptr, gctl::point3dc *op_ptr);
void gkernel_block_vz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, gctl::verbose_type_e verbose)
{
int i, j;
out_kernel.resize(obsp.size(), ele.size());
int o_size = obsp.size();
int e_size = ele.size();
gctl::progress_bar bar(o_size, "gkernel_vz");
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.at(i,j) = gkernel_block_vz_sig(ele.get(j), obsp.get(i));
}
}
return;
}
void gkernel_block_vzx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, gctl::verbose_type_e verbose)
{
int i, j;
out_kernel.resize(obsp.size(), ele.size());
int o_size = obsp.size();
int e_size = ele.size();
gctl::progress_bar bar(o_size, "gkernel_vzx");
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.at(i,j) = gkernel_block_vzx_sig(ele.get(j), obsp.get(i));
}
}
return;
}
void gkernel_block_vzy(gctl::matrix<double> &out_kernel, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, gctl::verbose_type_e verbose)
{
int i, j;
out_kernel.resize(obsp.size(), ele.size());
int o_size = obsp.size();
int e_size = ele.size();
gctl::progress_bar bar(o_size, "gkernel_vzy");
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.at(i,j) = gkernel_block_vzy_sig(ele.get(j), obsp.get(i));
}
}
return;
}
void gkernel_block_vzz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, gctl::verbose_type_e verbose)
{
int i, j;
out_kernel.resize(obsp.size(), ele.size());
int o_size = obsp.size();
int e_size = ele.size();
gctl::progress_bar bar(o_size, "gkernel_vzz");
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.at(i,j) = gkernel_block_vzz_sig(ele.get(j), obsp.get(i));
}
}
return;
}
void gobser_block_vz(gctl::array<double> &out_obs, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &rho,
gctl::verbose_type_e verbose)
{
int i, j;
out_obs.resize(obsp.size(), 0.0);
int o_size = obsp.size();
int e_size = ele.size();
gctl::progress_bar bar(e_size, "gobser_vz");
for (j = 0; j < e_size; j++)
{
if (verbose == gctl::FullMsg) bar.progressed(j);
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
if (rho[j] != 0.0)
{
#pragma omp parallel for private (i) schedule(guided)
for (i = 0; i < o_size; i++)
{
out_obs[i] += gkernel_block_vz_sig(ele.get(j), obsp.get(i)) * rho[j];
}
}
}
return;
}
void gobser_block_vzx(gctl::array<double> &out_obs, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &rho,
gctl::verbose_type_e verbose)
{
int i, j;
out_obs.resize(obsp.size(), 0.0);
int o_size = obsp.size();
int e_size = ele.size();
gctl::progress_bar bar(e_size, "gobser_vzx");
for (j = 0; j < e_size; j++)
{
if (verbose == gctl::FullMsg) bar.progressed(j);
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
if (rho[j] != 0.0)
{
#pragma omp parallel for private (i) schedule(guided)
for (i = 0; i < o_size; i++)
{
out_obs[i] += gkernel_block_vzx_sig(ele.get(j), obsp.get(i)) * rho[j];
}
}
}
return;
}
void gobser_block_vzy(gctl::array<double> &out_obs, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &rho,
gctl::verbose_type_e verbose)
{
int i, j;
out_obs.resize(obsp.size(), 0.0);
int o_size = obsp.size();
int e_size = ele.size();
gctl::progress_bar bar(e_size, "gobser_vzy");
for (j = 0; j < e_size; j++)
{
if (verbose == gctl::FullMsg) bar.progressed(j);
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
if (rho[j] != 0.0)
{
#pragma omp parallel for private (i) schedule(guided)
for (i = 0; i < o_size; i++)
{
out_obs[i] += gkernel_block_vzy_sig(ele.get(j), obsp.get(i)) * rho[j];
}
}
}
return;
}
void gobser_block_vzz(gctl::array<double> &out_obs, const gctl::array<gctl::block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &rho,
gctl::verbose_type_e verbose)
{
int i, j;
out_obs.resize(obsp.size(), 0.0);
int o_size = obsp.size();
int e_size = ele.size();
gctl::progress_bar bar(e_size, "gobser_vzz");
for (j = 0; j < e_size; j++)
{
if (verbose == gctl::FullMsg) bar.progressed(j);
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
if (rho[j] != 0.0)
{
#pragma omp parallel for private (i) schedule(guided)
for (i = 0; i < o_size; i++)
{
out_obs[i] += gkernel_block_vzz_sig(ele.get(j), obsp.get(i)) * rho[j];
}
}
}
return;
}
// define algorithm for individual element and observation point
double gkernel_block_vz_sig(gctl::block *ele_ptr, gctl::point3dc *op_ptr)
{
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
R222=sqrt((ele_ptr->ur->x - op_ptr->x)*(ele_ptr->ur->x - op_ptr->x)
+(ele_ptr->ur->y - op_ptr->y)*(ele_ptr->ur->y - op_ptr->y)
+(ele_ptr->ur->z - op_ptr->z)*(ele_ptr->ur->z - op_ptr->z));
R122=sqrt((ele_ptr->dl->x - op_ptr->x)*(ele_ptr->dl->x - op_ptr->x)
+(ele_ptr->ur->y - op_ptr->y)*(ele_ptr->ur->y - op_ptr->y)
+(ele_ptr->ur->z - op_ptr->z)*(ele_ptr->ur->z - op_ptr->z));
R212=sqrt((ele_ptr->ur->x - op_ptr->x)*(ele_ptr->ur->x - op_ptr->x)
+(ele_ptr->dl->y - op_ptr->y)*(ele_ptr->dl->y - op_ptr->y)
+(ele_ptr->ur->z - op_ptr->z)*(ele_ptr->ur->z - op_ptr->z));
R112=sqrt((ele_ptr->dl->x - op_ptr->x)*(ele_ptr->dl->x - op_ptr->x)
+(ele_ptr->dl->y - op_ptr->y)*(ele_ptr->dl->y - op_ptr->y)
+(ele_ptr->ur->z - op_ptr->z)*(ele_ptr->ur->z - op_ptr->z));
R221=sqrt((ele_ptr->ur->x - op_ptr->x)*(ele_ptr->ur->x - op_ptr->x)
+(ele_ptr->ur->y - op_ptr->y)*(ele_ptr->ur->y - op_ptr->y)
+(ele_ptr->dl->z - op_ptr->z)*(ele_ptr->dl->z - op_ptr->z));
R121=sqrt((ele_ptr->dl->x - op_ptr->x)*(ele_ptr->dl->x - op_ptr->x)
+(ele_ptr->ur->y - op_ptr->y)*(ele_ptr->ur->y - op_ptr->y)
+(ele_ptr->dl->z - op_ptr->z)*(ele_ptr->dl->z - op_ptr->z));
R211=sqrt((ele_ptr->ur->x - op_ptr->x)*(ele_ptr->ur->x - op_ptr->x)
+(ele_ptr->dl->y - op_ptr->y)*(ele_ptr->dl->y - op_ptr->y)
+(ele_ptr->dl->z - op_ptr->z)*(ele_ptr->dl->z - op_ptr->z));
R111=sqrt((ele_ptr->dl->x - op_ptr->x)*(ele_ptr->dl->x - op_ptr->x)
+(ele_ptr->dl->y - op_ptr->y)*(ele_ptr->dl->y - op_ptr->y)
+(ele_ptr->dl->z - op_ptr->z)*(ele_ptr->dl->z - op_ptr->z));
G222=(ele_ptr->ur->x-op_ptr->x)*log((ele_ptr->ur->y-op_ptr->y)+R222)
+(ele_ptr->ur->y-op_ptr->y)*log((ele_ptr->ur->x-op_ptr->x)+R222)
+(ele_ptr->ur->z-op_ptr->z)*gctl::arctg((ele_ptr->ur->z-op_ptr->z)
*R222/(ele_ptr->ur->x-op_ptr->x)/(ele_ptr->ur->y-op_ptr->y));
G122=(ele_ptr->dl->x-op_ptr->x)*log((ele_ptr->ur->y-op_ptr->y)+R122)
+(ele_ptr->ur->y-op_ptr->y)*log((ele_ptr->dl->x-op_ptr->x)+R122)
+(ele_ptr->ur->z-op_ptr->z)*gctl::arctg((ele_ptr->ur->z-op_ptr->z)
*R122/(ele_ptr->dl->x-op_ptr->x)/(ele_ptr->ur->y-op_ptr->y));
G212=(ele_ptr->ur->x-op_ptr->x)*log((ele_ptr->dl->y-op_ptr->y)+R212)
+(ele_ptr->dl->y-op_ptr->y)*log((ele_ptr->ur->x-op_ptr->x)+R212)
+(ele_ptr->ur->z-op_ptr->z)*gctl::arctg((ele_ptr->ur->z-op_ptr->z)
*R212/(ele_ptr->ur->x-op_ptr->x)/(ele_ptr->dl->y-op_ptr->y));
G112=(ele_ptr->dl->x-op_ptr->x)*log((ele_ptr->dl->y-op_ptr->y)+R112)
+(ele_ptr->dl->y-op_ptr->y)*log((ele_ptr->dl->x-op_ptr->x)+R112)
+(ele_ptr->ur->z-op_ptr->z)*gctl::arctg((ele_ptr->ur->z-op_ptr->z)
*R112/(ele_ptr->dl->x-op_ptr->x)/(ele_ptr->dl->y-op_ptr->y));
G221=(ele_ptr->ur->x-op_ptr->x)*log((ele_ptr->ur->y-op_ptr->y)+R221)
+(ele_ptr->ur->y-op_ptr->y)*log((ele_ptr->ur->x-op_ptr->x)+R221)
+(ele_ptr->dl->z-op_ptr->z)*gctl::arctg((ele_ptr->dl->z-op_ptr->z)
*R221/(ele_ptr->ur->x-op_ptr->x)/(ele_ptr->ur->y-op_ptr->y));
G121=(ele_ptr->dl->x-op_ptr->x)*log((ele_ptr->ur->y-op_ptr->y)+R121)
+(ele_ptr->ur->y-op_ptr->y)*log((ele_ptr->dl->x-op_ptr->x)+R121)
+(ele_ptr->dl->z-op_ptr->z)*gctl::arctg((ele_ptr->dl->z-op_ptr->z)
*R121/(ele_ptr->dl->x-op_ptr->x)/(ele_ptr->ur->y-op_ptr->y));
G211=(ele_ptr->ur->x-op_ptr->x)*log((ele_ptr->dl->y-op_ptr->y)+R211)
+(ele_ptr->dl->y-op_ptr->y)*log((ele_ptr->ur->x-op_ptr->x)+R211)
+(ele_ptr->dl->z-op_ptr->z)*gctl::arctg((ele_ptr->dl->z-op_ptr->z)
*R211/(ele_ptr->ur->x-op_ptr->x)/(ele_ptr->dl->y-op_ptr->y));
G111=(ele_ptr->dl->x-op_ptr->x)*log((ele_ptr->dl->y-op_ptr->y)+R111)
+(ele_ptr->dl->y-op_ptr->y)*log((ele_ptr->dl->x-op_ptr->x)+R111)
+(ele_ptr->dl->z-op_ptr->z)*gctl::arctg((ele_ptr->dl->z-op_ptr->z)
*R111/(ele_ptr->dl->x-op_ptr->x)/(ele_ptr->dl->y-op_ptr->y));
return 1.0e+8*GCTL_G0*(G222-G122-G212+G112-G221+G121+G211-G111);
}
double gkernel_block_vzx_sig(gctl::block *ele_ptr, gctl::point3dc *op_ptr)
{
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
R222=sqrt((ele_ptr->ur->x-op_ptr->x)*(ele_ptr->ur->x-op_ptr->x)
+(ele_ptr->ur->y-op_ptr->y)*(ele_ptr->ur->y-op_ptr->y)
+(ele_ptr->ur->z-op_ptr->z)*(ele_ptr->ur->z-op_ptr->z));
R122=sqrt((ele_ptr->dl->x-op_ptr->x)*(ele_ptr->dl->x-op_ptr->x)
+(ele_ptr->ur->y-op_ptr->y)*(ele_ptr->ur->y-op_ptr->y)
+(ele_ptr->ur->z-op_ptr->z)*(ele_ptr->ur->z-op_ptr->z));
R212=sqrt((ele_ptr->ur->x-op_ptr->x)*(ele_ptr->ur->x-op_ptr->x)
+(ele_ptr->dl->y-op_ptr->y)*(ele_ptr->dl->y-op_ptr->y)
+(ele_ptr->ur->z-op_ptr->z)*(ele_ptr->ur->z-op_ptr->z));
R112=sqrt((ele_ptr->dl->x-op_ptr->x)*(ele_ptr->dl->x-op_ptr->x)
+(ele_ptr->dl->y-op_ptr->y)*(ele_ptr->dl->y-op_ptr->y)
+(ele_ptr->ur->z-op_ptr->z)*(ele_ptr->ur->z-op_ptr->z));
R221=sqrt((ele_ptr->ur->x-op_ptr->x)*(ele_ptr->ur->x-op_ptr->x)
+(ele_ptr->ur->y-op_ptr->y)*(ele_ptr->ur->y-op_ptr->y)
+(ele_ptr->dl->z-op_ptr->z)*(ele_ptr->dl->z-op_ptr->z));
R121=sqrt((ele_ptr->dl->x-op_ptr->x)*(ele_ptr->dl->x-op_ptr->x)
+(ele_ptr->ur->y-op_ptr->y)*(ele_ptr->ur->y-op_ptr->y)
+(ele_ptr->dl->z-op_ptr->z)*(ele_ptr->dl->z-op_ptr->z));
R211=sqrt((ele_ptr->ur->x-op_ptr->x)*(ele_ptr->ur->x-op_ptr->x)
+(ele_ptr->dl->y-op_ptr->y)*(ele_ptr->dl->y-op_ptr->y)
+(ele_ptr->dl->z-op_ptr->z)*(ele_ptr->dl->z-op_ptr->z));
R111=sqrt((ele_ptr->dl->x-op_ptr->x)*(ele_ptr->dl->x-op_ptr->x)
+(ele_ptr->dl->y-op_ptr->y)*(ele_ptr->dl->y-op_ptr->y)
+(ele_ptr->dl->z-op_ptr->z)*(ele_ptr->dl->z-op_ptr->z));
G222=log((ele_ptr->ur->y-op_ptr->y)+R222);
G122=log((ele_ptr->ur->y-op_ptr->y)+R122);
G212=log((ele_ptr->dl->y-op_ptr->y)+R212);
G112=log((ele_ptr->dl->y-op_ptr->y)+R112);
G221=log((ele_ptr->ur->y-op_ptr->y)+R221);
G121=log((ele_ptr->ur->y-op_ptr->y)+R121);
G211=log((ele_ptr->dl->y-op_ptr->y)+R211);
G111=log((ele_ptr->dl->y-op_ptr->y)+R111);
return -1.0e+8*GCTL_G0*(G222-G122-G212+G112-G221+G121+G211-G111);
}
double gkernel_block_vzy_sig(gctl::block *ele_ptr, gctl::point3dc *op_ptr)
{
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
R222=sqrt((ele_ptr->ur->x-op_ptr->x)*(ele_ptr->ur->x-op_ptr->x)
+(ele_ptr->ur->y-op_ptr->y)*(ele_ptr->ur->y-op_ptr->y)
+(ele_ptr->ur->z-op_ptr->z)*(ele_ptr->ur->z-op_ptr->z));
R122=sqrt((ele_ptr->dl->x-op_ptr->x)*(ele_ptr->dl->x-op_ptr->x)
+(ele_ptr->ur->y-op_ptr->y)*(ele_ptr->ur->y-op_ptr->y)
+(ele_ptr->ur->z-op_ptr->z)*(ele_ptr->ur->z-op_ptr->z));
R212=sqrt((ele_ptr->ur->x-op_ptr->x)*(ele_ptr->ur->x-op_ptr->x)
+(ele_ptr->dl->y-op_ptr->y)*(ele_ptr->dl->y-op_ptr->y)
+(ele_ptr->ur->z-op_ptr->z)*(ele_ptr->ur->z-op_ptr->z));
R112=sqrt((ele_ptr->dl->x-op_ptr->x)*(ele_ptr->dl->x-op_ptr->x)
+(ele_ptr->dl->y-op_ptr->y)*(ele_ptr->dl->y-op_ptr->y)
+(ele_ptr->ur->z-op_ptr->z)*(ele_ptr->ur->z-op_ptr->z));
R221=sqrt((ele_ptr->ur->x-op_ptr->x)*(ele_ptr->ur->x-op_ptr->x)
+(ele_ptr->ur->y-op_ptr->y)*(ele_ptr->ur->y-op_ptr->y)
+(ele_ptr->dl->z-op_ptr->z)*(ele_ptr->dl->z-op_ptr->z));
R121=sqrt((ele_ptr->dl->x-op_ptr->x)*(ele_ptr->dl->x-op_ptr->x)
+(ele_ptr->ur->y-op_ptr->y)*(ele_ptr->ur->y-op_ptr->y)
+(ele_ptr->dl->z-op_ptr->z)*(ele_ptr->dl->z-op_ptr->z));
R211=sqrt((ele_ptr->ur->x-op_ptr->x)*(ele_ptr->ur->x-op_ptr->x)
+(ele_ptr->dl->y-op_ptr->y)*(ele_ptr->dl->y-op_ptr->y)
+(ele_ptr->dl->z-op_ptr->z)*(ele_ptr->dl->z-op_ptr->z));
R111=sqrt((ele_ptr->dl->x-op_ptr->x)*(ele_ptr->dl->x-op_ptr->x)
+(ele_ptr->dl->y-op_ptr->y)*(ele_ptr->dl->y-op_ptr->y)
+(ele_ptr->dl->z-op_ptr->z)*(ele_ptr->dl->z-op_ptr->z));
G222=log((ele_ptr->ur->x-op_ptr->x)+R222);
G122=log((ele_ptr->dl->x-op_ptr->x)+R122);
G212=log((ele_ptr->ur->x-op_ptr->x)+R212);
G112=log((ele_ptr->dl->x-op_ptr->x)+R112);
G221=log((ele_ptr->ur->x-op_ptr->x)+R221);
G121=log((ele_ptr->dl->x-op_ptr->x)+R121);
G211=log((ele_ptr->ur->x-op_ptr->x)+R211);
G111=log((ele_ptr->dl->x-op_ptr->x)+R111);
return -1.0e+8*GCTL_G0*(G222-G122-G212+G112-G221+G121+G211-G111);
}
double gkernel_block_vzz_sig(gctl::block *ele_ptr, gctl::point3dc *op_ptr)
{
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
R222=sqrt((ele_ptr->ur->x-op_ptr->x)*(ele_ptr->ur->x-op_ptr->x)
+(ele_ptr->ur->y-op_ptr->y)*(ele_ptr->ur->y-op_ptr->y)
+(ele_ptr->ur->z-op_ptr->z)*(ele_ptr->ur->z-op_ptr->z));
R122=sqrt((ele_ptr->dl->x-op_ptr->x)*(ele_ptr->dl->x-op_ptr->x)
+(ele_ptr->ur->y-op_ptr->y)*(ele_ptr->ur->y-op_ptr->y)
+(ele_ptr->ur->z-op_ptr->z)*(ele_ptr->ur->z-op_ptr->z));
R212=sqrt((ele_ptr->ur->x-op_ptr->x)*(ele_ptr->ur->x-op_ptr->x)
+(ele_ptr->dl->y-op_ptr->y)*(ele_ptr->dl->y-op_ptr->y)
+(ele_ptr->ur->z-op_ptr->z)*(ele_ptr->ur->z-op_ptr->z));
R112=sqrt((ele_ptr->dl->x-op_ptr->x)*(ele_ptr->dl->x-op_ptr->x)
+(ele_ptr->dl->y-op_ptr->y)*(ele_ptr->dl->y-op_ptr->y)
+(ele_ptr->ur->z-op_ptr->z)*(ele_ptr->ur->z-op_ptr->z));
R221=sqrt((ele_ptr->ur->x-op_ptr->x)*(ele_ptr->ur->x-op_ptr->x)
+(ele_ptr->ur->y-op_ptr->y)*(ele_ptr->ur->y-op_ptr->y)
+(ele_ptr->dl->z-op_ptr->z)*(ele_ptr->dl->z-op_ptr->z));
R121=sqrt((ele_ptr->dl->x-op_ptr->x)*(ele_ptr->dl->x-op_ptr->x)
+(ele_ptr->ur->y-op_ptr->y)*(ele_ptr->ur->y-op_ptr->y)
+(ele_ptr->dl->z-op_ptr->z)*(ele_ptr->dl->z-op_ptr->z));
R211=sqrt((ele_ptr->ur->x-op_ptr->x)*(ele_ptr->ur->x-op_ptr->x)
+(ele_ptr->dl->y-op_ptr->y)*(ele_ptr->dl->y-op_ptr->y)
+(ele_ptr->dl->z-op_ptr->z)*(ele_ptr->dl->z-op_ptr->z));
R111=sqrt((ele_ptr->dl->x-op_ptr->x)*(ele_ptr->dl->x-op_ptr->x)
+(ele_ptr->dl->y-op_ptr->y)*(ele_ptr->dl->y-op_ptr->y)
+(ele_ptr->dl->z-op_ptr->z)*(ele_ptr->dl->z-op_ptr->z));
G222=atan((ele_ptr->ur->x-op_ptr->x)*(ele_ptr->ur->y-op_ptr->y)
/(R222*(ele_ptr->ur->z-op_ptr->z)));
G122=atan((ele_ptr->dl->x-op_ptr->x)*(ele_ptr->ur->y-op_ptr->y)
/(R122*(ele_ptr->ur->z-op_ptr->z)));
G212=atan((ele_ptr->ur->x-op_ptr->x)*(ele_ptr->dl->y-op_ptr->y)
/(R212*(ele_ptr->ur->z-op_ptr->z)));
G112=atan((ele_ptr->dl->x-op_ptr->x)*(ele_ptr->dl->y-op_ptr->y)
/(R112*(ele_ptr->ur->z-op_ptr->z)));
G221=atan((ele_ptr->ur->x-op_ptr->x)*(ele_ptr->ur->y-op_ptr->y)
/(R221*(ele_ptr->dl->z-op_ptr->z)));
G121=atan((ele_ptr->dl->x-op_ptr->x)*(ele_ptr->ur->y-op_ptr->y)
/(R121*(ele_ptr->dl->z-op_ptr->z)));
G211=atan((ele_ptr->ur->x-op_ptr->x)*(ele_ptr->dl->y-op_ptr->y)
/(R211*(ele_ptr->dl->z-op_ptr->z)));
G111=atan((ele_ptr->dl->x-op_ptr->x)*(ele_ptr->dl->y-op_ptr->y)
/(R111*(ele_ptr->dl->z-op_ptr->z)));
return -1.0e+8*GCTL_G0*(G222-G122-G212+G112-G221+G121+G211-G111);
}
Reference in New Issue View Git Blame Copy Permalink