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

update source

Browse Source
This commit is contained in:
张壹
2021-07-18 22:12:09 +08:00
parent c3efc2928d
commit 266169698e
7 changed files with 410 additions and 410 deletions
Download Patch File Download Diff File Expand all files Collapse all files

38
src/magtet.cpp
View File

@@ -17,7 +17,7 @@ void magtet::read_tet(std::string tet_name)
void magtet::read_magz(std::string mag_name)
{
std::vector<gctl::point3d_c> tmp_vec;
std::vector<gctl::point3dc> tmp_vec;
gctl::read_text2vector(mag_name, "", tmp_vec);
if (tmp_vec.size() != ele_num_)
@@ -35,7 +35,7 @@ void magtet::init_site(std::string para)
// try to use the para as a file name
if (access(para.c_str(), F_OK) != -1)
{
std::vector<gctl::point3d_c> tmp_vec;
std::vector<gctl::point3dc> tmp_vec;
gctl::read_text2vector(para, "", tmp_vec);
site_num_ = tmp_vec.size();
@@ -138,7 +138,7 @@ void magtet::cal_tensors()
enorm_.resize(ele_num_, 12);
etang_.resize(ele_num_, 12);
gctl::point3d_c v1, v2, v3, nf, ne;
gctl::point3dc v1, v2, v3, nf, ne;
for (int e = 0; e < ele_num_; ++e)
{
for (int i = 0; i < 4; ++i)
@@ -160,12 +160,12 @@ void magtet::cal_tensors()
return;
}
double magtet::mag_potential(const gctl::tetrahedron &tet, const gctl::point3d_c &site,
const gctl::point3d_c &mz, gctl::point3d_c *fn, gctl::point3d_c *en, gctl::point3d_c *et)
double magtet::mag_potential(const gctl::tetrahedron &tet, const gctl::point3dc &site,
const gctl::point3dc &mz, gctl::point3dc *fn, gctl::point3dc *en, gctl::point3dc *et)
{
double Rij_minus, Rij_plus, Sij_plus, Sij_minus, Rij0, mij0, wi0;
double part1, part2, k0, absw, mdotn, beta;
gctl::point3d_c oi;
gctl::point3dc oi;
double out_pot = 0.0;
for (int f = 0; f < 4; ++f)
@@ -212,14 +212,14 @@ double magtet::mag_potential(const gctl::tetrahedron &tet, const gctl::point3d_c
return out_pot;
}
gctl::point3d_c magtet::mag_gradient(const gctl::tetrahedron &tet, const gctl::point3d_c &site,
const gctl::point3d_c &mz, gctl::point3d_c *fn, gctl::point3d_c *en, gctl::point3d_c *et)
gctl::point3dc magtet::mag_gradient(const gctl::tetrahedron &tet, const gctl::point3dc &site,
const gctl::point3dc &mz, gctl::point3dc *fn, gctl::point3dc *en, gctl::point3dc *et)
{
double Rij_minus, Rij_plus, Sij_plus, Sij_minus, Rij0, mij0, wi0;
double mdotn, beta, Aij, sig, absw;
gctl::point3d_c oi, k1, part1, part2;
gctl::point3dc oi, k1, part1, part2;
gctl::point3d_c out_grad(0.0, 0.0, 0.0);
gctl::point3dc out_grad(0.0, 0.0, 0.0);
for (int f = 0; f < 4; ++f)
{
k1.set(0.0, 0.0, 0.0);
@@ -278,17 +278,17 @@ gctl::point3d_c magtet::mag_gradient(const gctl::tetrahedron &tet, const gctl::p
return out_grad;
}
gctl::tensor magtet::mag_tensor(const gctl::tetrahedron &tet, const gctl::point3d_c &site,
const gctl::point3d_c &mz, gctl::point3d_c *fn, gctl::point3d_c *en, gctl::point3d_c *et)
gctl::tensor magtet::mag_tensor(const gctl::tetrahedron &tet, const gctl::point3dc &site,
const gctl::point3dc &mz, gctl::point3dc *fn, gctl::point3dc *en, gctl::point3dc *et)
{
double Rij_minus, Rij_plus, Sij_plus, Sij_minus, Rij0, mij0, wi0;
double mdotn, beta, sig, absw;
double factor_n_mij, factor_tij;
gctl::point3d_c oi, grad_Aij;
gctl::point3d_c grad_Rij_plus, grad_Rij_minus, grad_Sij_plus, grad_Sij_minus;
gctl::point3d_c grad_mij0, grad_Rij0, grad_abs_wi0;
gctl::point3d_c grad_a_plus, grad_b_plus, grad_a_minus, grad_b_minus;
gctl::point3d_c grad_betaij_plus, grad_betaij_minus, grad_betaij;
gctl::point3dc oi, grad_Aij;
gctl::point3dc grad_Rij_plus, grad_Rij_minus, grad_Sij_plus, grad_Sij_minus;
gctl::point3dc grad_mij0, grad_Rij0, grad_abs_wi0;
gctl::point3dc grad_a_plus, grad_b_plus, grad_a_minus, grad_b_minus;
gctl::point3dc grad_betaij_plus, grad_betaij_minus, grad_betaij;
double a_plus, b_plus, a_minus, b_minus;
double k3;
gctl::tensor tmp_k, k2;
@@ -404,7 +404,7 @@ void magtet::total_gradient()
{
int e, s;
// malloc space
mag_grad_.resize(site_num_, gctl::point3d_c(0.0, 0.0, 0.0));
mag_grad_.resize(site_num_, gctl::point3dc(0.0, 0.0, 0.0));
gctl::utility::progress_bar bar(ele_num_, "MagGradient");
for (e = 0; e < ele_num_; ++e)
@@ -491,7 +491,7 @@ void magtet::routine(const char *para_file)
if (!matched)
{
throw gctl::runtime_error("Invalid calculating type. From magtet::main(...)");
throw gctl::runtime_error("Invalid calculating type. From magtet::routine(...)");
}
// output results

26
src/magtet.h
View File

@@ -31,12 +31,12 @@ public:
void write_text(std::string out_name);
void cal_tensors();
double mag_potential(const gctl::tetrahedron &tet, const gctl::point3d_c &site, const gctl::point3d_c &mz,
gctl::point3d_c *fn, gctl::point3d_c *en, gctl::point3d_c *et);
gctl::point3d_c mag_gradient(const gctl::tetrahedron &tet, const gctl::point3d_c &site, const gctl::point3d_c &mz,
gctl::point3d_c *fn, gctl::point3d_c *en, gctl::point3d_c *et);
gctl::tensor mag_tensor(const gctl::tetrahedron &tet, const gctl::point3d_c &site, const gctl::point3d_c &mz,
gctl::point3d_c *fn, gctl::point3d_c *en, gctl::point3d_c *et);
double mag_potential(const gctl::tetrahedron &tet, const gctl::point3dc &site, const gctl::point3dc &mz,
gctl::point3dc *fn, gctl::point3dc *en, gctl::point3dc *et);
gctl::point3dc mag_gradient(const gctl::tetrahedron &tet, const gctl::point3dc &site, const gctl::point3dc &mz,
gctl::point3dc *fn, gctl::point3dc *en, gctl::point3dc *et);
gctl::tensor mag_tensor(const gctl::tetrahedron &tet, const gctl::point3dc &site, const gctl::point3dc &mz,
gctl::point3dc *fn, gctl::point3dc *en, gctl::point3dc *et);
void total_potential();
void total_gradient();
void total_tensor();
@@ -47,16 +47,16 @@ public:
protected:
int node_num_, ele_num_, site_num_;
gctl::array<gctl::vertex3d_c> node_; // 四面体元素的顶点集
gctl::array<gctl::vertex3dc> node_; // 四面体元素的顶点集
gctl::array<gctl::tetrahedron> ele_; // 四面体的元素集
gctl::array<gctl::point3d_c> magz_; // 磁化矢量
gctl::array2d<gctl::point3d_c> fnorm_; // 面的外法线矢量
gctl::array2d<gctl::point3d_c> enorm_; // 边的外法线矢量
gctl::array2d<gctl::point3d_c> etang_;
gctl::array<gctl::point3d_c> site_; // 观测点位置集
gctl::array<gctl::point3dc> magz_; // 磁化矢量
gctl::array2d<gctl::point3dc> fnorm_; // 面的外法线矢量
gctl::array2d<gctl::point3dc> enorm_; // 边的外法线矢量
gctl::array2d<gctl::point3dc> etang_;
gctl::array<gctl::point3dc> site_; // 观测点位置集
gctl::array<double> mag_pot_; // 正演磁位
gctl::array<gctl::point3d_c> mag_grad_; // 正演磁场分量
gctl::array<gctl::point3dc> mag_grad_; // 正演磁场分量
gctl::array<gctl::tensor> mag_tensor_; // 正演磁场张量
gctl::utility::getoption gopt_; // 参数文件解释对象