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

update code

Browse Source
This commit is contained in:
张壹 2021-08-06 21:37:28 +08:00
parent 393ffdee91
commit f3e58b8fd7
8 changed files with 6388 additions and 6383 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
data/prism/prism_B.nc
View File

Binary file not shown.

BIN
data/prism/prism_T.nc
View File

Binary file not shown.

BIN
data/prism/prism_V.nc
View File

Binary file not shown.

7492
data/prism/prism_mag_grad.txt
View File

File diff suppressed because it is too large Load Diff

5218
data/prism/prism_mag_pot.txt
View File

File diff suppressed because it is too large Load Diff

8
log.txt
View File

@ -1,5 +1,5 @@
tet-file = data/pipeline/pipeline tet-file = data/prism/prism.1
mag-file = data/pipeline/magz.txt mag-file = (0,0,200)
site-file = data/pipeline/site.txt site-file = -30/30/-30/30/10/81/81
obs-file = data/pipeline/pipeline obs-file = data/prism/prism
cal-type = potential gradient tensor cal-type = potential gradient tensor

31
src/magtet.cpp
View File

@ -1,10 +1,5 @@
#include "magtet.h" #include "magtet.h"
magtet::magtet(){}
magtet::~magtet(){}
void magtet::read_tet(std::string tet_name) void magtet::read_tet(std::string tet_name)
{ {
gctl::read_Tetgen_node(tet_name, node_); gctl::read_Tetgen_node(tet_name, node_);
@ -130,44 +125,46 @@ void magtet::write_text(std::string out_name)
void magtet::cal_tensors() void magtet::cal_tensors()
{ {
ele_para_.resize(ele_num_); //ele_para_.resize(ele_num_);
gctl::point3dc v1, v2, v3, nf, ne; gctl::point3dc v1, v2, v3, nf, ne;
for (int e = 0; e < ele_num_; ++e) for (int e = 0; e < ele_num_; ++e)
{ {
ele_[e].att = new magtet_para;
for (int i = 0; i < 4; ++i) for (int i = 0; i < 4; ++i)
{ {
v1 = *ele_[e].fget(i, 1) - *ele_[e].fget(i, 0); v1 = *ele_[e].fget(i, 1) - *ele_[e].fget(i, 0);
v2 = *ele_[e].fget(i, 2) - *ele_[e].fget(i, 0); v2 = *ele_[e].fget(i, 2) - *ele_[e].fget(i, 0);
nf = gctl::cross(v1, v2).normal(); nf = gctl::cross(v1, v2).normal();
// The space is declared by the read_magz() function // The space is declared by the read_magz() function
ele_para_[e].fnorm[i] = nf; ele_[e].att->fnorm[i] = nf;
for (int j = 0; j < 3; ++j) for (int j = 0; j < 3; ++j)
{ {
v3 = *ele_[e].fget(i, (j+1)%3) - *ele_[e].fget(i, j); v3 = *ele_[e].fget(i, (j+1)%3) - *ele_[e].fget(i, j);
ne = gctl::cross(v3, nf).normal(); ne = gctl::cross(v3, nf).normal();
ele_para_[e].enorm[j+i*3] = ne; ele_[e].att->enorm[j+i*3] = ne;
ele_para_[e].etang[j+i*3] = gctl::cross(nf, ne); ele_[e].att->etang[j+i*3] = gctl::cross(nf, ne);
} }
ele_para_[e].mag_amp[i] = gctl::dot(magz_[e], ele_para_[e].fnorm[i]); ele_[e].att->mag_amp[i] = gctl::dot(magz_[e], ele_[e].att->fnorm[i]);
} }
// link magtet_para to tetrahedron's attribute // link magtet_para to tetrahedron's attribute
ele_[e].att = ele_para_.get(e); //ele_[e].att = ele_para_.get(e);
} }
return; return;
} }
double magtet::mag_potential(const gctl::tetrahedron &tet, const gctl::point3dc &site) double magtet::mag_potential(const mag_tetrahedron &tet, const gctl::point3dc &site)
{ {
double Rij_minus, Rij_plus, Sij_plus, Sij_minus, Rij0, mij0, wi0; double Rij_minus, Rij_plus, Sij_plus, Sij_minus, Rij0, mij0, wi0;
double part1, part2, k0, absw, beta; double part1, part2, k0, absw, beta;
gctl::point3dc oi; gctl::point3dc oi;
// get attribute pointer // get attribute pointer
magtet_para *mpara = (magtet_para*)tet.att; magtet_para *mpara = tet.att;
double out_pot = 0.0; double out_pot = 0.0;
for (int f = 0; f < 4; ++f) for (int f = 0; f < 4; ++f)
@ -213,14 +210,14 @@ double magtet::mag_potential(const gctl::tetrahedron &tet, const gctl::point3dc
return out_pot; return out_pot;
} }
gctl::point3dc magtet::mag_gradient(const gctl::tetrahedron &tet, const gctl::point3dc &site) gctl::point3dc magtet::mag_gradient(const mag_tetrahedron &tet, const gctl::point3dc &site)
{ {
double Rij_minus, Rij_plus, Sij_plus, Sij_minus, Rij0, mij0, wi0; double Rij_minus, Rij_plus, Sij_plus, Sij_minus, Rij0, mij0, wi0;
double beta, Aij, sig, absw; double beta, Aij, sig, absw;
gctl::point3dc oi, k1, part1, part2; gctl::point3dc oi, k1, part1, part2;
// get attribute pointer // get attribute pointer
magtet_para *mpara = (magtet_para*)tet.att; magtet_para *mpara = tet.att;
gctl::point3dc 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) for (int f = 0; f < 4; ++f)
@ -280,7 +277,7 @@ gctl::point3dc magtet::mag_gradient(const gctl::tetrahedron &tet, const gctl::po
return out_grad; return out_grad;
} }
gctl::tensor magtet::mag_tensor(const gctl::tetrahedron &tet, const gctl::point3dc &site) gctl::tensor magtet::mag_tensor(const mag_tetrahedron &tet, const gctl::point3dc &site)
{ {
double Rij_minus, Rij_plus, Sij_plus, Sij_minus, Rij0, mij0, wi0; double Rij_minus, Rij_plus, Sij_plus, Sij_minus, Rij0, mij0, wi0;
double beta, sig, absw; double beta, sig, absw;
@ -295,7 +292,7 @@ gctl::tensor magtet::mag_tensor(const gctl::tetrahedron &tet, const gctl::point3
gctl::tensor tmp_k, k2; gctl::tensor tmp_k, k2;
// get attribute pointer // get attribute pointer
magtet_para *mpara = (magtet_para*)tet.att; magtet_para *mpara = tet.att;
gctl::tensor out_tensor(0.0); gctl::tensor out_tensor(0.0);
for (int f = 0; f < 4; ++f) for (int f = 0; f < 4; ++f)

22
src/magtet.h
View File

@ -25,11 +25,19 @@ struct magtet_para
gctl::point3dc fnorm[4], enorm[12], etang[12]; gctl::point3dc fnorm[4], enorm[12], etang[12];
}; };
typedef gctl::type_tetrahedron<magtet_para> mag_tetrahedron;
class magtet class magtet
{ {
public: public:
magtet(); magtet(){}
virtual ~magtet(); virtual ~magtet()
{
for (int i = 0; i < ele_num_; ++i)
{
delete ele_[i].att;
}
}
void read_tet(std::string tet_name); void read_tet(std::string tet_name);
void init_magz(std::string para); void init_magz(std::string para);
@ -37,9 +45,9 @@ public:
void write_text(std::string out_name); void write_text(std::string out_name);
void cal_tensors(); void cal_tensors();
double mag_potential(const gctl::tetrahedron &tet, const gctl::point3dc &site); double mag_potential(const mag_tetrahedron &tet, const gctl::point3dc &site);
gctl::point3dc mag_gradient(const gctl::tetrahedron &tet, const gctl::point3dc &site); gctl::point3dc mag_gradient(const mag_tetrahedron &tet, const gctl::point3dc &site);
gctl::tensor mag_tensor(const gctl::tetrahedron &tet, const gctl::point3dc &site); gctl::tensor mag_tensor(const mag_tetrahedron &tet, const gctl::point3dc &site);
void total_potential(); void total_potential();
void total_gradient(); void total_gradient();
void total_tensor(); void total_tensor();
@ -50,9 +58,9 @@ public:
protected: protected:
int node_num_, ele_num_, site_num_; int node_num_, ele_num_, site_num_;
gctl::array<gctl::vertex3dc> node_; // 四面体元素的顶点集 gctl::array<gctl::vertex3dc> node_; // 四面体元素的顶点集
gctl::array<gctl::tetrahedron> ele_; // 四面体的元素集 gctl::array<mag_tetrahedron> ele_; // 四面体的元素集
gctl::array<gctl::point3dc> magz_; // 磁化矢量数组 gctl::array<gctl::point3dc> magz_; // 磁化矢量数组
gctl::array<magtet_para> ele_para_; // 四面体的张量属性与磁化矢量 //gctl::array<magtet_para> ele_para_; // 四面体的张量属性与磁化矢量
gctl::array<gctl::point3dc> site_; // 观测点位置集 gctl::array<gctl::point3dc> site_; // 观测点位置集
gctl::array<double> mag_pot_; // 正演磁位 gctl::array<double> mag_pot_; // 正演磁位