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-22 10:31:29 +08:00
parent d50a0b0c17
commit 44e1ab4d51
13 changed files with 13394 additions and 13397 deletions
Download Patch File Download Diff File Expand all files Collapse all files

30
src/magtet.cpp
View File

@@ -28,12 +28,7 @@ void magtet::init_magz(std::string para)
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];
}
magz_.resize(tmp_vec);
gctl::destroy_vector(tmp_vec);
return;
}
@@ -42,10 +37,10 @@ void magtet::init_magz(std::string para)
gctl::point3dc tmp_magz;
tmp_magz.str(para);
ele_para_.resize(ele_num_);
magz_.resize(ele_num_);
for (int i = 0; i < ele_num_; ++i)
{
ele_para_[i].magz = tmp_magz;
magz_[i] = tmp_magz;
}
return;
}
@@ -135,6 +130,8 @@ void magtet::write_text(std::string out_name)
void magtet::cal_tensors()
{
ele_para_.resize(ele_num_);
gctl::point3dc v1, v2, v3, nf, ne;
for (int e = 0; e < ele_num_; ++e)
{
@@ -153,6 +150,8 @@ void magtet::cal_tensors()
ele_para_[e].enorm[j+i*3] = ne;
ele_para_[e].etang[j+i*3] = gctl::cross(nf, ne);
}
ele_para_[e].mag_amp[i] = gctl::dot(magz_[e], ele_para_[e].fnorm[i]);
}
// link magtet_para to tetrahedron's attribute
@@ -164,7 +163,7 @@ void magtet::cal_tensors()
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;
double part1, part2, k0, absw, beta;
gctl::point3dc oi;
// get attribute pointer
@@ -208,8 +207,7 @@ double magtet::mag_potential(const gctl::tetrahedron &tet, const gctl::point3dc
k0 += (part1 - part2);
}
mdotn = gctl::dot(mpara->magz, mpara->fnorm[f]);
out_pot += Cm*k0*mdotn;
out_pot += Cm*k0*mpara->mag_amp[f];
}
return out_pot;
@@ -218,7 +216,7 @@ double magtet::mag_potential(const gctl::tetrahedron &tet, const gctl::point3dc
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;
double beta, Aij, sig, absw;
gctl::point3dc oi, k1, part1, part2;
// get attribute pointer
@@ -276,8 +274,7 @@ gctl::point3dc magtet::mag_gradient(const gctl::tetrahedron &tet, const gctl::po
k1 = k1 - (part1 + part2);
}
mdotn = gctl::dot(mpara->magz, mpara->fnorm[f]);
out_grad = out_grad - 1e+9*Cm*mdotn*k1;
out_grad = out_grad - 1e+9*Cm*k1*mpara->mag_amp[f];
}
return out_grad;
@@ -286,7 +283,7 @@ gctl::point3dc magtet::mag_gradient(const gctl::tetrahedron &tet, const gctl::po
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;
double beta, sig, absw;
double factor_n_mij, factor_tij;
gctl::point3dc oi, grad_Aij;
gctl::point3dc grad_Rij_plus, grad_Rij_minus, grad_Sij_plus, grad_Sij_minus;
@@ -380,8 +377,7 @@ gctl::tensor magtet::mag_tensor(const gctl::tetrahedron &tet, const gctl::point3
k2 = k2 - k3*tmp_k;
}
mdotn = gctl::dot(mpara->magz, mpara->fnorm[f]);
out_tensor = out_tensor - 1e+9*Cm*mdotn*k2;
out_tensor = out_tensor - 1e+9*Cm*k2*mpara->mag_amp[f];
}
return out_tensor;

3
src/magtet.h
View File

@@ -21,7 +21,7 @@
struct magtet_para
{
gctl::point3dc magz;
double mag_amp[4];
gctl::point3dc fnorm[4], enorm[12], etang[12];
};
@@ -51,6 +51,7 @@ protected:
int node_num_, ele_num_, site_num_;
gctl::array<gctl::vertex3dc> node_; // 四面体元素的顶点集
gctl::array<gctl::tetrahedron> ele_; // 四面体的元素集
gctl::array<gctl::point3dc> magz_; // 磁化矢量数组
gctl::array<magtet_para> ele_para_; // 四面体的张量属性与磁化矢量
gctl::array<gctl::point3dc> site_; // 观测点位置集