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

update src

Browse Source
This commit is contained in:
张壹
2021-07-21 20:25:30 +08:00
parent 47d623e37d
commit 2b832d8830
18 changed files with 1676 additions and 970 deletions
Download Patch File Download Diff File Expand all files Collapse all files

143
src/magtet.cpp
View File

@@ -15,18 +15,37 @@ void magtet::read_tet(std::string tet_name)
return;
}
void magtet::read_magz(std::string mag_name)
void magtet::init_magz(std::string para)
{
std::vector<gctl::point3dc> tmp_vec;
gctl::read_text2vector(mag_name, "", tmp_vec);
if (tmp_vec.size() != ele_num_)
// try to use the para as a file name
if (access(para.c_str(), F_OK) != -1)
{
throw gctl::length_error("Invalid magnetization vector size. From magtet::read_magz(...)");
std::vector<gctl::point3dc> tmp_vec;
gctl::read_text2vector(para, "", tmp_vec);
if (tmp_vec.size() != ele_num_)
{
throw gctl::length_error("Invalid magnetization vector size. From magtet::init_magz(...)");
}
ele_para_.resize(tmp_vec.size());
for (int i = 0; i < tmp_vec.size(); ++i)
{
ele_para_[i].magz = tmp_vec[i];
}
gctl::destroy_vector(tmp_vec);
}
magz_.resize(tmp_vec);
gctl::destroy_vector(tmp_vec);
// try to use the para in the format <mag_x>,<mag_y>,<mag_z>
gctl::point3dc tmp_magz;
tmp_magz.str(para);
ele_para_.resize(ele_num_);
for (int i = 0; i < ele_num_; ++i)
{
ele_para_[i].magz = tmp_magz;
}
return;
}
@@ -115,11 +134,6 @@ void magtet::write_text(std::string out_name)
void magtet::cal_tensors()
{
// malloc space
fnorm_.resize(ele_num_, 4);
enorm_.resize(ele_num_, 12);
etang_.resize(ele_num_, 12);
gctl::point3dc v1, v2, v3, nf, ne;
for (int e = 0; e < ele_num_; ++e)
{
@@ -128,27 +142,33 @@ void magtet::cal_tensors()
v1 = *ele_[e].fget(i, 1) - *ele_[e].fget(i, 0);
v2 = *ele_[e].fget(i, 2) - *ele_[e].fget(i, 0);
nf = gctl::cross(v1, v2).normal();
fnorm_[e][i] = nf;
// The space is declared by the read_magz() function
ele_para_[e].fnorm[i] = nf;
for (int j = 0; j < 3; ++j)
{
v3 = *ele_[e].fget(i, (j+1)%3) - *ele_[e].fget(i, j);
ne = gctl::cross(v3, nf).normal();
enorm_[e][j+i*3] = ne;
etang_[e][j+i*3] = gctl::cross(nf, ne);
ele_para_[e].enorm[j+i*3] = ne;
ele_para_[e].etang[j+i*3] = gctl::cross(nf, ne);
}
}
// link magtet_para to tetrahedron's attribute
ele_[e].att = ele_para_.get(e);
}
return;
}
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 magtet::mag_potential(const gctl::tetrahedron &tet, const gctl::point3dc &site)
{
double Rij_minus, Rij_plus, Sij_plus, Sij_minus, Rij0, mij0, wi0;
double part1, part2, k0, absw, mdotn, beta;
gctl::point3dc oi;
// get attribute pointer
magtet_para *mpara = (magtet_para*)tet.att;
double out_pot = 0.0;
for (int f = 0; f < 4; ++f)
{
@@ -160,14 +180,14 @@ double magtet::mag_potential(const gctl::tetrahedron &tet, const gctl::point3dc
if (j == 0)
{
wi0 = gctl::dot(site - *tet.fget(f, j), fn[f]);
wi0 = gctl::dot(site - *tet.fget(f, j), mpara->fnorm[f]);
absw = std::abs(wi0);
}
oi = site - wi0*fn[f];
Sij_minus = gctl::dot(*tet.fget(f, j) - oi, et[3*f+j]);
Sij_plus = gctl::dot(*tet.fget(f, (j+1)%3) - oi, et[3*f+j]);
mij0 = gctl::dot(*tet.fget(f, j) - oi, en[3*f+j]);
oi = site - wi0*mpara->fnorm[f];
Sij_minus = gctl::dot(*tet.fget(f, j) - oi, mpara->etang[3*f+j]);
Sij_plus = gctl::dot(*tet.fget(f, (j+1)%3) - oi, mpara->etang[3*f+j]);
mij0 = gctl::dot(*tet.fget(f, j) - oi, mpara->enorm[3*f+j]);
Rij0 = std::sqrt(wi0*wi0 + mij0*mij0);
part2 = 0.0;
@@ -187,20 +207,22 @@ double magtet::mag_potential(const gctl::tetrahedron &tet, const gctl::point3dc
k0 += (part1 - part2);
}
mdotn = gctl::dot(mz, fn[f]);
mdotn = gctl::dot(mpara->magz, mpara->fnorm[f]);
out_pot += Cm*k0*mdotn;
}
return out_pot;
}
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)
gctl::point3dc magtet::mag_gradient(const gctl::tetrahedron &tet, const gctl::point3dc &site)
{
double Rij_minus, Rij_plus, Sij_plus, Sij_minus, Rij0, mij0, wi0;
double mdotn, beta, Aij, sig, absw;
gctl::point3dc oi, k1, part1, part2;
// get attribute pointer
magtet_para *mpara = (magtet_para*)tet.att;
gctl::point3dc out_grad(0.0, 0.0, 0.0);
for (int f = 0; f < 4; ++f)
{
@@ -212,15 +234,15 @@ gctl::point3dc magtet::mag_gradient(const gctl::tetrahedron &tet, const gctl::po
if (j == 0)
{
wi0 = gctl::dot(site - *tet.fget(f, j), fn[f]);
wi0 = gctl::dot(site - *tet.fget(f, j), mpara->fnorm[f]);
sig = sign(wi0);
absw = std::abs(wi0);
}
oi = site - wi0*fn[f];
Sij_minus = gctl::dot(*tet.fget(f, j) - oi, et[3*f+j]);
Sij_plus = gctl::dot(*tet.fget(f, (j+1)%3) - oi, et[3*f+j]);
mij0 = gctl::dot(*tet.fget(f, j) - oi, en[3*f+j]);
oi = site - wi0*mpara->fnorm[f];
Sij_minus = gctl::dot(*tet.fget(f, j) - oi, mpara->etang[3*f+j]);
Sij_plus = gctl::dot(*tet.fget(f, (j+1)%3) - oi, mpara->etang[3*f+j]);
mij0 = gctl::dot(*tet.fget(f, j) - oi, mpara->enorm[3*f+j]);
Rij0 = std::sqrt(wi0*wi0 + mij0*mij0);
part2.set(0.0, 0.0, 0.0);
@@ -229,7 +251,7 @@ gctl::point3dc magtet::mag_gradient(const gctl::tetrahedron &tet, const gctl::po
beta = atan((mij0*Sij_plus)/(Rij0*Rij0 + absw*Rij_plus))
- atan((mij0*Sij_minus)/(Rij0*Rij0 + absw*Rij_minus));
part2 = sig*beta*fn[f];
part2 = sig*beta*mpara->fnorm[f];
}
if (std::abs(Rij0) > TOL)
@@ -249,19 +271,18 @@ gctl::point3dc magtet::mag_gradient(const gctl::tetrahedron &tet, const gctl::po
throw gctl::runtime_error("Observation site on edge. From magtet::mag_gradient()");
}
part1 = Aij * en[3*f+j];
part1 = Aij * mpara->enorm[3*f+j];
k1 = k1 - (part1 + part2);
}
mdotn = gctl::dot(mz, fn[f]);
mdotn = gctl::dot(mpara->magz, mpara->fnorm[f]);
out_grad = out_grad - 1e+9*Cm*mdotn*k1;
}
return out_grad;
}
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)
gctl::tensor magtet::mag_tensor(const gctl::tetrahedron &tet, const gctl::point3dc &site)
{
double Rij_minus, Rij_plus, Sij_plus, Sij_minus, Rij0, mij0, wi0;
double mdotn, beta, sig, absw;
@@ -275,6 +296,9 @@ gctl::tensor magtet::mag_tensor(const gctl::tetrahedron &tet, const gctl::point3
double k3;
gctl::tensor tmp_k, k2;
// get attribute pointer
magtet_para *mpara = (magtet_para*)tet.att;
gctl::tensor out_tensor(0.0);
for (int f = 0; f < 4; ++f)
{
@@ -287,43 +311,43 @@ gctl::tensor magtet::mag_tensor(const gctl::tetrahedron &tet, const gctl::point3
if (j == 0)
{
wi0 = gctl::dot(site - *tet.fget(f, j), fn[f]);
wi0 = gctl::dot(site - *tet.fget(f, j), mpara->fnorm[f]);
sig = sign(wi0);
absw = std::abs(wi0);
}
oi = site - wi0*fn[f];
Sij_minus = gctl::dot(*tet.fget(f, j) - oi, et[3*f+j]);
Sij_plus = gctl::dot(*tet.fget(f, (j+1)%3) - oi, et[3*f+j]);
mij0 = gctl::dot(*tet.fget(f, j) - oi, en[3*f+j]);
oi = site - wi0*mpara->fnorm[f];
Sij_minus = gctl::dot(*tet.fget(f, j) - oi, mpara->etang[3*f+j]);
Sij_plus = gctl::dot(*tet.fget(f, (j+1)%3) - oi, mpara->etang[3*f+j]);
mij0 = gctl::dot(*tet.fget(f, j) - oi, mpara->enorm[3*f+j]);
Rij0 = std::sqrt(wi0*wi0 + mij0*mij0);
if (std::abs(Rij0) > TOL)
{
factor_n_mij = -1.0*(Sij_plus/(Rij0*Rij0*Rij_plus) - Sij_minus/(Rij0*Rij0*Rij_minus));
factor_tij = -1.0/Rij_plus + 1.0/Rij_minus;
grad_Aij = (wi0*factor_n_mij)*fn[f] + factor_tij*et[3*f+j]
- (mij0*factor_n_mij)*en[3*f+j];
grad_Aij = (wi0*factor_n_mij)*mpara->fnorm[f] + factor_tij*mpara->etang[3*f+j]
- (mij0*factor_n_mij)*mpara->enorm[3*f+j];
}
else
{
factor_tij = -1.0/Rij_plus + 1.0/Rij_minus;
grad_Aij = factor_tij*et[3*f+j];
grad_Aij = factor_tij*mpara->etang[3*f+j];
}
//tmp_k = gctl::kron(grad_Aij, enorm_[e][3*f+j]);
tmp_k = gctl::kron(en[3*f+j], grad_Aij);
tmp_k = gctl::kron(mpara->enorm[3*f+j], grad_Aij);
k2 = k2 - tmp_k;
if (absw > TOL)
{
grad_Rij_plus = (1.0/Rij_plus)*(site - *tet.fget(f, (j+1)%3));
grad_Rij_minus = (1.0/Rij_minus)*(site - *tet.fget(f, j));
grad_Sij_plus = -1.0*et[3*f+j];
grad_Sij_minus = -1.0*et[3*f+j];
grad_mij0 = -1.0*en[3*f+j];
grad_Rij0 = (1.0/Rij0)*(wi0*fn[f] - mij0*en[3*f+j]);
grad_abs_wi0 = sig*fn[f];
grad_Sij_plus = -1.0*mpara->etang[3*f+j];
grad_Sij_minus = -1.0*mpara->etang[3*f+j];
grad_mij0 = -1.0*mpara->enorm[3*f+j];
grad_Rij0 = (1.0/Rij0)*(wi0*mpara->fnorm[f] - mij0*mpara->enorm[3*f+j]);
grad_abs_wi0 = sig*mpara->fnorm[f];
a_plus = Rij0*Rij0 + absw*Rij_plus;
b_plus = mij0*Sij_plus;
grad_a_plus = (2.0*Rij0)*grad_Rij0 + Rij_plus*grad_abs_wi0 + absw*grad_Rij_plus;
@@ -337,7 +361,7 @@ gctl::tensor magtet::mag_tensor(const gctl::tetrahedron &tet, const gctl::point3
grad_betaij = grad_betaij_plus - grad_betaij_minus;
tmp_k = gctl::kron(grad_betaij, fn[f]);
tmp_k = gctl::kron(grad_betaij, mpara->fnorm[f]);
k2 = k2 - sig*tmp_k;
}
else
@@ -351,11 +375,11 @@ gctl::tensor magtet::mag_tensor(const gctl::tetrahedron &tet, const gctl::point3
if (k3 != 0.0)
{
tmp_k = gctl::kron(fn[f], fn[f]);
tmp_k = gctl::kron(mpara->fnorm[f], mpara->fnorm[f]);
k2 = k2 - k3*tmp_k;
}
mdotn = gctl::dot(mz, fn[f]);
mdotn = gctl::dot(mpara->magz, mpara->fnorm[f]);
out_tensor = out_tensor - 1e+9*Cm*mdotn*k2;
}
@@ -376,7 +400,7 @@ void magtet::total_potential()
#pragma omp parallel for private (s) schedule (guided)
for (s = 0; s < site_num_; ++s)
{
mag_pot_[s] += mag_potential(ele_[e], site_[s], magz_[e], fnorm_.get(e), enorm_.get(e), etang_.get(e));
mag_pot_[s] += mag_potential(ele_[e], site_[s]);
}
}
return;
@@ -396,8 +420,7 @@ void magtet::total_gradient()
#pragma omp parallel for private (s) schedule (guided)
for (s = 0; s < site_num_; ++s)
{
mag_grad_[s] = mag_grad_[s] +
mag_gradient(ele_[e], site_[s], magz_[e], fnorm_.get(e), enorm_.get(e), etang_.get(e));
mag_grad_[s] = mag_grad_[s] + mag_gradient(ele_[e], site_[s]);
}
}
return;
@@ -417,8 +440,7 @@ void magtet::total_tensor()
#pragma omp parallel for private (s) schedule (guided)
for (s = 0; s < site_num_; ++s)
{
mag_tensor_[s] = mag_tensor_[s] +
mag_tensor(ele_[e], site_[s], magz_[e], fnorm_.get(e), enorm_.get(e), etang_.get(e));
mag_tensor_[s] = mag_tensor_[s] + mag_tensor(ele_[e], site_[s]);
}
}
return;
@@ -440,7 +462,7 @@ void magtet::routine(const char *para_file)
// read tetgen files
RUN_ECHO(read_tet(gopt_.get_value(TETFILE)), "Reading 3D model file");
// read magnetization file
RUN_ECHO(read_magz(gopt_.get_value(MAGFILE)), "Reading magnetization file");
RUN_ECHO(init_magz(gopt_.get_value(MAGFILE)), "Reading magnetization file");
// read site file
RUN_ECHO(init_site(gopt_.get_value(SITEFILE)), "Initiating observations points");
// initialize tensors
@@ -492,7 +514,8 @@ int main(int argc, char const *argv[])
std::clog << "<option> = <value>\n";
std::clog << "Available options:\n";
std::clog << "tet-file: 3D model file. this option needs both .node and .ele files generated by Tetgen. No file extension is needed\n";
std::clog << "mag-file: magnetization vectors. Line format: <magz-x> <magz-y> <magz-z>\n";
std::clog << "mag-file: magnetization vectors. format1: <filename> format2: (<magz-x>,<magz-y>,<magz-z>)\n";
std::clog << "mag-file: line format: <magz-x> <magz-y> <magz-z>\n";
std::clog << "site-file: observation points. format1: <filename> format2: <xmin>/<xmax>/<ymin>/<ymax>/<z>/<xnum>/<ynum>\n";
std::clog << "site-file: line format: <x> <y> <z>\n";
std::clog << "obs-file: output observation files. No file extension is needed\n";

25
src/magtet.h
View File

@@ -17,7 +17,13 @@
#define OBSFILE "obs-file"
#define CALTYPE "cal-type" // potential, gradient or tensor
#define RUN_ECHO(action, msg) do {std::clog << msg << " ... \n"; action;} while(0);
#define RUN_ECHO(action, msg) do {std::clog << msg << "... \n"; action;} while(0);
struct magtet_para
{
gctl::point3dc magz;
gctl::point3dc fnorm[4], enorm[12], etang[12];
};
class magtet
{
@@ -26,33 +32,26 @@ public:
virtual ~magtet();
void read_tet(std::string tet_name);
void read_magz(std::string mag_name);
void init_magz(std::string para);
void init_site(std::string para);
void write_text(std::string out_name);
void cal_tensors();
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);
double mag_potential(const gctl::tetrahedron &tet, const gctl::point3dc &site);
gctl::point3dc mag_gradient(const gctl::tetrahedron &tet, const gctl::point3dc &site);
gctl::tensor mag_tensor(const gctl::tetrahedron &tet, const gctl::point3dc &site);
void total_potential();
void total_gradient();
void total_tensor();
double sign(double s);
void routine(const char *para_file);
protected:
int node_num_, ele_num_, site_num_;
gctl::array<gctl::vertex3dc> node_; // 四面体元素的顶点集
gctl::array<gctl::tetrahedron> ele_; // 四面体的元素集
gctl::array<gctl::point3dc> magz_; // 磁化矢量
gctl::array2d<gctl::point3dc> fnorm_; // 面的外法线矢量
gctl::array2d<gctl::point3dc> enorm_; // 边的外法线矢量
gctl::array2d<gctl::point3dc> etang_;
gctl::array<magtet_para> ele_para_; // 四面体的张量属性与磁化矢量
gctl::array<gctl::point3dc> site_; // 观测点位置集
gctl::array<double> mag_pot_; // 正演磁位