zhangyiss/gm3d-gui
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

add function for forward mdling using meshtool3d mesh and phy files

Browse Source
This commit is contained in:
kiciollo 2019-09-12 18:55:37 +08:00
parent ab2a765c91
commit 716939c93a
31 changed files with 289 additions and 227 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/CMakeLists.txt
View File

@ -13,6 +13,8 @@ add_executable(gm3d_gui WIN32 MACOSX_BUNDLE ${DIR_SRC})
include_directories(/usr/local/include)
find_library(FLTK_LIBRARY fltk HINTS /usr/local/lib)
find_library(HINTS /usr/local/lib)
#find_library(FLTK_LIBRARY fltk HINTS /usr/local/lib)
target_link_libraries(gm3d_gui PUBLIC ${FLTK_LIBRARY})

8
src/add_interface_block.cpp
View File

@ -1,5 +1,5 @@
#include "gm3d.h"
//我们读入一个界面数据 插值计算每个块体中心位置的值 然后按情况赋值
//我们读入一个界面数据 插值计算每个块体中心位置的值 然后按情况赋值
int GM3D::AddInterfaceBlock(modelist para_list){
int m, n, xnum, ynum;
double xs,xe,xmin,xmax,dx;
@ -18,7 +18,7 @@ int GM3D::AddInterfaceBlock(modelist para_list){
if (open_infile(infile,filename)) return -1;
while (getline(infile,temp_str)){
//#range必须出现在数据之前
//#range必须出现在数据之前
if (*(temp_str.begin()) == '#'){
if (6 == sscanf(temp_str.c_str(),"# range=%lf/%lf/%lf/%lf/%lf/%lf",&xs,&dx,&xe,&ys,&dy,&ye)){
xmin = MIN(xs,xe); xmax = MAX(xs,xe);
@ -75,7 +75,7 @@ int GM3D::AddInterfaceBlock(modelist para_list){
if (model_cube_[i].cen.z >= temp_topo.z){
if (model_block_val_[i] == BDL_MAX)
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
else
model_block_val_[i] += para_list.mod_value;
model_added = true;
@ -92,7 +92,7 @@ int GM3D::AddInterfaceBlock(modelist para_list){
if (model_cube_[i].cen.z <= temp_topo.z){
if (model_block_val_[i] == BDL_MAX)
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
else
model_block_val_[i] += para_list.mod_value;
model_added = true;

8
src/add_model_gui.cpp
View File

@ -18,7 +18,7 @@ Fl_Return_Button *can_add_btn=(Fl_Return_Button *)0;
Fl_Button *sig_add_btn=(Fl_Button *)0;
Fl_Input *mod_val_input=(Fl_Input *)0;
// 声明三个字符串用于保存新建的模型类型和赋值类型
// 声明三个字符串用于保存新建的模型类型和赋值类型
char add_mod_type[1024];
char add_val_type[1024];
char agn_val_type[1024];
@ -81,8 +81,8 @@ void cb_bot_val_rbtn(Fl_Button*, void*){
return;
}
// 声明一个字符串用于保存新建的模型参数并添加至列表中
// 模型参数的排列为<模型类型> <赋值类型> <物理参数字符串> <几何参数字符串>
// 声明一个字符串用于保存新建的模型参数并添加至列表中
// 模型参数的排列为<模型类型> <赋值类型> <物理参数字符串> <几何参数字符串>
char add_mod_para[1024];
void cb_sig_add_btn(Fl_Widget*, void*){
@ -101,7 +101,7 @@ void cb_sig_add_btn(Fl_Widget*, void*){
mod_para_brw->add(add_mod_para);
return;
}
//取消添加模型
//取消添加模型
void cb_can_add_btn(Fl_Widget*, void*){
add_mod_win->hide();
return;

2
src/add_models.cpp
View File

@ -10,7 +10,7 @@ int GM3D::AddModels(char* filename){
while(getline(infile,temp_str)){
if (*(temp_str.begin()) == '#') continue;
else{
//按每行5个数据解析 初始化为含观测值与不确定度的观测点
//按每行5个数据解析 初始化为含观测值与不确定度的观测点
if (4 == sscanf(temp_str.c_str(),"%s %s %lf %s",
temp_list.mod_type,temp_list.val_type,&temp_list.mod_value,temp_list.mod_para)){
model_list_.push_back(temp_list);

6
src/add_regular_block.cpp
View File

@ -16,7 +16,7 @@ int GM3D::AddRegularBlock(modelist para_list){
if (model_cube_[i].cen.x >= xmin && model_cube_[i].cen.x <= xmax &&
model_cube_[i].cen.y >= ymin && model_cube_[i].cen.y <= ymax &&
model_cube_[i].cen.z >= zmin && model_cube_[i].cen.z <= zmax){
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
model_added = true;
}
}
@ -27,7 +27,7 @@ int GM3D::AddRegularBlock(modelist para_list){
model_cube_[i].cen.y >= ymin && model_cube_[i].cen.y <= ymax &&
model_cube_[i].cen.z >= zmin && model_cube_[i].cen.z <= zmax){
if (model_block_val_[i] == BDL_MAX)
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
else
model_block_val_[i] += para_list.mod_value;
model_added = true;
@ -39,7 +39,7 @@ int GM3D::AddRegularBlock(modelist para_list){
if (model_cube_[i].cen.x >= xmin && model_cube_[i].cen.x <= xmax &&
model_cube_[i].cen.y >= ymin && model_cube_[i].cen.y <= ymax &&
model_cube_[i].cen.z >= zmin && model_cube_[i].cen.z <= zmax){
model_block_val_[i] = BDL_MAX; //注意重复赋值的块体会覆盖
model_block_val_[i] = BDL_MAX; //注意重复赋值的块体会覆盖
model_added = true;
}
}

6
src/add_sphere_block.cpp
View File

@ -20,7 +20,7 @@ int GM3D::AddSphereBlock(modelist para_list){
rad_limit = rad_x*rad_y*rad_z/sqrt(pow(rad_y*rad_z*sin(theta)*cos(phi),2) + pow(rad_x*rad_z*sin(theta)*sin(phi),2) + pow(rad_x*rad_y*cos(theta),2));
if (dist <= rad_limit){
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
model_added = true;
}
}
@ -39,7 +39,7 @@ int GM3D::AddSphereBlock(modelist para_list){
if (dist <= rad_limit){
if (model_block_val_[i] == BDL_MAX)
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
else
model_block_val_[i] += para_list.mod_value;
model_added = true;
@ -59,7 +59,7 @@ int GM3D::AddSphereBlock(modelist para_list){
rad_limit = rad_x*rad_y*rad_z/sqrt(pow(rad_y*rad_z*sin(theta)*cos(phi),2) + pow(rad_x*rad_z*sin(theta)*sin(phi),2) + pow(rad_x*rad_y*cos(theta),2));
if (dist <= rad_limit){
model_block_val_[i] = BDL_MAX; //注意重复赋值的块体会覆盖
model_block_val_[i] = BDL_MAX; //注意重复赋值的块体会覆盖
model_added = true;
}
}

12
src/add_tilted_block.cpp
View File

@ -19,7 +19,7 @@ int GM3D::AddTiltedBlock(modelist para_list){
if (!strcmp(para_list.val_type,"replace")){
for (int i = 0; i < model_num_; i++){
//计算当前层的x与y范围
//计算当前层的x与y范围
layer_xmin = (model_cube_[i].cen.z - zmin)*(xmin_2 - xmin_1)/(zmax - zmin) + xmin_1;
layer_xmax = (model_cube_[i].cen.z - zmin)*(xmax_2 - xmax_1)/(zmax - zmin) + xmax_1;
layer_ymin = (model_cube_[i].cen.z - zmin)*(ymin_2 - ymin_1)/(zmax - zmin) + ymin_1;
@ -28,14 +28,14 @@ int GM3D::AddTiltedBlock(modelist para_list){
if (model_cube_[i].cen.x >= layer_xmin && model_cube_[i].cen.x <= layer_xmax &&
model_cube_[i].cen.y >= layer_ymin && model_cube_[i].cen.y <= layer_ymax &&
model_cube_[i].cen.z >= zmin && model_cube_[i].cen.z <= zmax){
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
model_added = true;
}
}
}
else if (!strcmp(para_list.val_type,"add")){
for (int i = 0; i < model_num_; i++){
//计算当前层的x与y范围
//计算当前层的x与y范围
layer_xmin = (model_cube_[i].cen.z - zmin)*(xmin_2 - xmin_1)/(zmax - zmin) + xmin_1;
layer_xmax = (model_cube_[i].cen.z - zmin)*(xmax_2 - xmax_1)/(zmax - zmin) + xmax_1;
layer_ymin = (model_cube_[i].cen.z - zmin)*(ymin_2 - ymin_1)/(zmax - zmin) + ymin_1;
@ -45,7 +45,7 @@ int GM3D::AddTiltedBlock(modelist para_list){
model_cube_[i].cen.y >= layer_ymin && model_cube_[i].cen.y <= layer_ymax &&
model_cube_[i].cen.z >= zmin && model_cube_[i].cen.z <= zmax){
if (model_block_val_[i] == BDL_MAX)
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
else
model_block_val_[i] += para_list.mod_value;
model_added = true;
@ -54,7 +54,7 @@ int GM3D::AddTiltedBlock(modelist para_list){
}
else if (!strcmp(para_list.val_type,"erase")){
for (int i = 0; i < model_num_; i++){
//计算当前层的x与y范围
//计算当前层的x与y范围
layer_xmin = (model_cube_[i].cen.z - zmin)*(xmin_2 - xmin_1)/(zmax - zmin) + xmin_1;
layer_xmax = (model_cube_[i].cen.z - zmin)*(xmax_2 - xmax_1)/(zmax - zmin) + xmax_1;
layer_ymin = (model_cube_[i].cen.z - zmin)*(ymin_2 - ymin_1)/(zmax - zmin) + ymin_1;
@ -63,7 +63,7 @@ int GM3D::AddTiltedBlock(modelist para_list){
if (model_cube_[i].cen.x >= layer_xmin && model_cube_[i].cen.x <= layer_xmax &&
model_cube_[i].cen.y >= layer_ymin && model_cube_[i].cen.y <= layer_ymax &&
model_cube_[i].cen.z >= zmin && model_cube_[i].cen.z <= zmax){
model_block_val_[i] = BDL_MAX; //注意重复赋值的块体会覆盖
model_block_val_[i] = BDL_MAX; //注意重复赋值的块体会覆盖
model_added = true;
}
}

20
src/build_regular_grid.cpp
View File

@ -46,26 +46,26 @@ int GM3D::BuildRegularGrid(char* space_para){
while(x >= xmin && x <= xmax){
z = zs;
while(z >= zmin && z <= zmax){
//添加denMod
//添加denMod
temp_cu.cen.id = model_cube_.size();
temp_cu.cen.x = x; temp_cu.cen.y = y; temp_cu.cen.z = z;
//添加mshVert
//添加mshVert
for (int i = 0; i < 8; i++){
temp_cp.id = model_vert_.size(); //添加msh的顶点索引为mshVert的大小
temp_cp.x = temp_cu.cen.x - sign[i][0]*temp_cu.dx; //左下底角
temp_cp.id = model_vert_.size(); //添加msh的顶点索引为mshVert的大小
temp_cp.x = temp_cu.cen.x - sign[i][0]*temp_cu.dx; //左下底角
temp_cp.y = temp_cu.cen.y - sign[i][1]*temp_cu.dy;
temp_cp.z = temp_cu.cen.z - sign[i][2]*temp_cu.dz;
temp_id_str = cpoint_id(temp_cp);
imsp = map_str_point.find(temp_id_str);
//利用map_vert查到当前顶点是否存在,这里需要注意，如果顶点已经存在则只需要将顶点索引置为已存在顶点的索引，不增加顶点计数
//利用map_vert查到当前顶点是否存在,这里需要注意，如果顶点已经存在则只需要将顶点索引置为已存在顶点的索引，不增加顶点计数
if(imsp!=map_str_point.end()){
temp_cu.ids[i] = imsp->second.id;
}
//若为新的顶点则将其增加到两个映射和一个链表中
//若为新的顶点则将其增加到两个映射和一个链表中
else{
temp_cu.ids[i] = temp_cp.id;//新的顶点索引等于顶点集的数量
model_vert_.push_back(temp_cp);//将新产生的顶点保存到顶点链表中
map_str_point[temp_id_str] = temp_cp;//将新产生的顶点保存到顶点位置映射中
temp_cu.ids[i] = temp_cp.id;//新的顶点索引等于顶点集的数量
model_vert_.push_back(temp_cp);//将新产生的顶点保存到顶点链表中
map_str_point[temp_id_str] = temp_cp;//将新产生的顶点保存到顶点位置映射中
}
}
model_cube_.push_back(temp_cu);
@ -82,7 +82,7 @@ int GM3D::BuildRegularGrid(char* space_para){
else{
vert_num_ = model_vert_.size();
model_num_ = model_cube_.size();
model_block_val_.resize(model_num_,BDL_MAX); //初始化模型块体值为BDL_MAX
model_block_val_.resize(model_num_,BDL_MAX); //初始化模型块体值为BDL_MAX
}
return 0;
}

8
src/forward_delta_t.cpp
View File

@ -10,18 +10,18 @@ int GM3D::ForwardDeltaT(char* noise_level,char* mag_para){
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
//初始化正演值和不确定度
//初始化正演值和不确定度
for (int i = 0; i < obs_num_; i++){
obs_p_[i].val = obs_p_[i].dev = 0.0;
}
//确定噪声水平
//确定噪声水平
double noise_mean, noise_dev, temp_noise;
if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
noise_mean = noise_dev = 0.0;
}
//确定磁化参数
//确定磁化参数
if (4 != sscanf(mag_para,"%lf/%lf/%lf/%lf",&I0,&A0,&I,&A)){
I0 = I = 90; A0 = A = 0;
}
@ -39,7 +39,7 @@ int GM3D::ForwardDeltaT(char* noise_level,char* mag_para){
k6=-sin(I0)*sin(I);
}
//添加高斯噪声值
//添加高斯噪声值
default_random_engine generator;
normal_distribution<double> dist(noise_mean, noise_dev);

8
src/forward_delta_tx.cpp
View File

@ -10,18 +10,18 @@ int GM3D::ForwardDeltaTx(char* noise_level,char* mag_para){
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
//初始化正演值和不确定度
//初始化正演值和不确定度
for (int i = 0; i < obs_num_; i++){
obs_p_[i].val = obs_p_[i].dev = 0.0;
}
//确定噪声水平
//确定噪声水平
double noise_mean, noise_dev, temp_noise;
if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
noise_mean = noise_dev = 0.0;
}
//确定磁化参数
//确定磁化参数
if (4 != sscanf(mag_para,"%lf/%lf/%lf/%lf",&I0,&A0,&I,&A)){
I0 = I = 90; A0 = A = 0;
}
@ -39,7 +39,7 @@ int GM3D::ForwardDeltaTx(char* noise_level,char* mag_para){
k6=-sin(I0)*sin(I);
}
//添加高斯噪声值
//添加高斯噪声值
default_random_engine generator;
normal_distribution<double> dist(noise_mean, noise_dev);

8
src/forward_delta_ty.cpp
View File

@ -10,18 +10,18 @@ int GM3D::ForwardDeltaTy(char* noise_level,char* mag_para){
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
//初始化正演值和不确定度
//初始化正演值和不确定度
for (int i = 0; i < obs_num_; i++){
obs_p_[i].val = obs_p_[i].dev = 0.0;
}
//确定噪声水平
//确定噪声水平
double noise_mean, noise_dev, temp_noise;
if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
noise_mean = noise_dev = 0.0;
}
//确定磁化参数
//确定磁化参数
if (4 != sscanf(mag_para,"%lf/%lf/%lf/%lf",&I0,&A0,&I,&A)){
I0 = I = 90; A0 = A = 0;
}
@ -39,7 +39,7 @@ int GM3D::ForwardDeltaTy(char* noise_level,char* mag_para){
k6=-sin(I0)*sin(I);
}
//添加高斯噪声值
//添加高斯噪声值
default_random_engine generator;
normal_distribution<double> dist(noise_mean, noise_dev);

8
src/forward_delta_tz.cpp
View File

@ -11,18 +11,18 @@ int GM3D::ForwardDeltaTz(char* noise_level,char* mag_para){
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
//初始化正演值和不确定度
//初始化正演值和不确定度
for (int i = 0; i < obs_num_; i++){
obs_p_[i].val = obs_p_[i].dev = 0.0;
}
//确定噪声水平
//确定噪声水平
double noise_mean, noise_dev, temp_noise;
if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
noise_mean = noise_dev = 0.0;
}
//确定磁化参数
//确定磁化参数
if (4 != sscanf(mag_para,"%lf/%lf/%lf/%lf",&I0,&A0,&I,&A)){
I0 = I = 90; A0 = A = 0;
}
@ -40,7 +40,7 @@ int GM3D::ForwardDeltaTz(char* noise_level,char* mag_para){
k6=-sin(I0)*sin(I);
}
//添加高斯噪声值
//添加高斯噪声值
default_random_engine generator;
normal_distribution<double> dist(noise_mean, noise_dev);

8
src/forward_hax.cpp
View File

@ -10,18 +10,18 @@ int GM3D::ForwardHax(char* noise_level,char* mag_para){
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
//初始化正演值和不确定度
//初始化正演值和不确定度
for (int i = 0; i < obs_num_; i++){
obs_p_[i].val = obs_p_[i].dev = 0.0;
}
//确定噪声水平
//确定噪声水平
double noise_mean, noise_dev, temp_noise;
if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
noise_mean = noise_dev = 0.0;
}
//确定磁化参数
//确定磁化参数
if (4 != sscanf(mag_para,"%lf/%lf/%lf/%lf",&I0,&A0,&I,&A)){
I = 90; A = 0;
}
@ -34,7 +34,7 @@ int GM3D::ForwardHax(char* noise_level,char* mag_para){
Gamma=sin(I);
}
//添加高斯噪声值
//添加高斯噪声值
default_random_engine generator;
normal_distribution<double> dist(noise_mean, noise_dev);

8
src/forward_hay.cpp
View File

@ -10,18 +10,18 @@ int GM3D::ForwardHay(char* noise_level,char* mag_para){
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
//初始化正演值和不确定度
//初始化正演值和不确定度
for (int i = 0; i < obs_num_; i++){
obs_p_[i].val = obs_p_[i].dev = 0.0;
}
//确定噪声水平
//确定噪声水平
double noise_mean, noise_dev, temp_noise;
if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
noise_mean = noise_dev = 0.0;
}
//确定磁化参数
//确定磁化参数
if (4 != sscanf(mag_para,"%lf/%lf/%lf/%lf",&I0,&A0,&I,&A)){
I = 90; A = 0;
}
@ -34,7 +34,7 @@ int GM3D::ForwardHay(char* noise_level,char* mag_para){
Gamma=sin(I);
}
//添加高斯噪声值
//添加高斯噪声值
default_random_engine generator;
normal_distribution<double> dist(noise_mean, noise_dev);

6
src/forward_vz.cpp
View File

@ -6,18 +6,18 @@ int GM3D::ForwardVz(char* noise_level){
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
//初始化正演值和不确定度
//初始化正演值和不确定度
for (int i = 0; i < obs_num_; i++){
obs_p_[i].val = obs_p_[i].dev = 0.0;
}
//确定噪声水平
//确定噪声水平
double noise_mean, noise_dev, temp_noise;
if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
noise_mean = noise_dev = 0.0;
}
//添加高斯噪声值
//添加高斯噪声值
default_random_engine generator;
normal_distribution<double> dist(noise_mean, noise_dev);

4
src/forward_vzx.cpp
View File

@ -6,7 +6,7 @@ int GM3D::ForwardVzx(char* noise_level){
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
//初始化正演值和不确定度
//初始化正演值和不确定度
for (int i = 0; i < obs_num_; i++){
obs_p_[i].val = obs_p_[i].dev = 0.0;
}
@ -16,7 +16,7 @@ int GM3D::ForwardVzx(char* noise_level){
noise_mean = noise_dev = 0.0;
}
//添加高斯噪声值
//添加高斯噪声值
default_random_engine generator;
normal_distribution<double> dist(noise_mean, noise_dev);

4
src/forward_vzy.cpp
View File

@ -6,7 +6,7 @@ int GM3D::ForwardVzy(char* noise_level){
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
//初始化正演值和不确定度
//初始化正演值和不确定度
for (int i = 0; i < obs_num_; i++){
obs_p_[i].val = obs_p_[i].dev = 0.0;
}
@ -16,7 +16,7 @@ int GM3D::ForwardVzy(char* noise_level){
noise_mean = noise_dev = 0.0;
}
//添加高斯噪声值
//添加高斯噪声值
default_random_engine generator;
normal_distribution<double> dist(noise_mean, noise_dev);

4
src/forward_vzz.cpp
View File

@ -6,7 +6,7 @@ int GM3D::ForwardVzz(char* noise_level){
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
//初始化正演值和不确定度
//初始化正演值和不确定度
for (int i = 0; i < obs_num_; i++){
obs_p_[i].val = obs_p_[i].dev = 0.0;
}
@ -16,7 +16,7 @@ int GM3D::ForwardVzz(char* noise_level){
noise_mean = noise_dev = 0.0;
}
//添加高斯噪声值
//添加高斯噪声值
default_random_engine generator;
normal_distribution<double> dist(noise_mean, noise_dev);

8
src/forward_za.cpp
View File

@ -10,18 +10,18 @@ int GM3D::ForwardZa(char* noise_level,char* mag_para){
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
//初始化正演值和不确定度
//初始化正演值和不确定度
for (int i = 0; i < obs_num_; i++){
obs_p_[i].val = obs_p_[i].dev = 0.0;
}
//确定噪声水平
//确定噪声水平
double noise_mean, noise_dev, temp_noise;
if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
noise_mean = noise_dev = 0.0;
}
//确定磁化参数
//确定磁化参数
if (4 != sscanf(mag_para,"%lf/%lf/%lf/%lf",&I0,&A0,&I,&A)){
I = 90; A = 0;
}
@ -34,7 +34,7 @@ int GM3D::ForwardZa(char* noise_level,char* mag_para){
Gamma=sin(I);
}
//添加高斯噪声值
//添加高斯噪声值
default_random_engine generator;
normal_distribution<double> dist(noise_mean, noise_dev);

39
src/gm3d.h
View File

@ -7,23 +7,24 @@ class GM3D{
public:
GM3D(){}
~GM3D(){}
int BuildRegularGrid(char*); //初始化反演模型空间
int AddModels(char*); //读取模型块体参数文件
void AddModels_GUI(); //添加模型块体参数
int AddRegularBlock(modelist); //添加普通模型块体
int AddTiltedBlock(modelist); //添加倾斜模型块体
int AddSphereBlock(modelist); //添加球体椭球体块体
int AddInterfaceBlock(modelist); //添加密度界面
//模型操作
int BuildRegularGrid(char*); //初始化反演模型空间
int AddModels(char*); //读取模型块体参数文件
void AddModels_GUI(); //添加模型块体参数
int AddRegularBlock(modelist); //添加普通模型块体
int AddTiltedBlock(modelist); //添加倾斜模型块体
int AddSphereBlock(modelist); //添加球体椭球体块体
int AddInterfaceBlock(modelist); //添加密度界面
//模型操作
int ReadModel(char*,char*);
//输出模型
int RegisteredOuput(bool); //注册输出的块体模型
int OutMshFile(char*,string); //输出模型文件
int OutNeighborFile(char*,char*); //输出模型块体或顶点的相邻关系 暂缓
//观测数据
int ReadModel_mst(char*);
//输出模型
int RegisteredOuput(bool); //注册输出的块体模型
int OutMshFile(char*,string); //输出模型文件
int OutNeighborFile(char*,char*); //输出模型块体或顶点的相邻关系 暂缓
//观测数据
int InitObs(char*);
int OutObs(char*);
//正演函数
//正演函数
int ForwardVz(char*);
int ForwardVzx(char*);
int ForwardVzy(char*);
@ -35,16 +36,16 @@ public:
int ForwardHax(char*,char*);
int ForwardHay(char*,char*);
int ForwardZa(char*,char*);
//gm3d_gui函数
//gm3d_gui函数
void get_model_list(modelistArray);
private:
int obs_num_, model_num_, vert_num_;
//正演数组
//正演数组
obspointArray obs_p_;
_2dArray input_models_;
_1sArray input_model_names_;
_1dArray forward_model_;
//模型数据
//模型数据
cubeArray model_cube_;
cpointArray model_vert_;
_1dArray model_block_val_;
@ -55,7 +56,7 @@ private:
_int2intMap vert_out_map_;
_int2intMap ele_data_out_map_;
_2iArray model_vert_neighbor_; //暂缓
_2iArray model_cube_neighbor_; //暂缓
_2iArray model_vert_neighbor_; //暂缓
_2iArray model_cube_neighbor_; //暂缓
};
#endif

48
src/gm3d_gui.cpp
View File

@ -64,7 +64,7 @@ void cb_mesh_file_btn(Fl_Widget*, void*){
// Block until user picks something.
while(chooser.shown()) { Fl::wait(); }
//从文件的绝对路径中分离文件名并从输出框输出
//从文件的绝对路径中分离文件名并从输出框输出
if (chooser.value() != NULL)
{
strcpy(mesh_filename,chooser.value());
@ -92,7 +92,7 @@ void cb_mod_para_file_btn(Fl_Widget*, void*){
// Block until user picks something.
while(chooser.shown()) { Fl::wait(); }
//下面我们读入文件并添加模型参数到列表显示
//下面我们读入文件并添加模型参数到列表显示
std::ifstream modpara_in;
modpara_in.open(chooser.value());
if (!modpara_in)
@ -149,7 +149,7 @@ void cb_mod_file_out_btn(Fl_Widget*, void*){
// Block until user picks something.
while(chooser.shown()) { Fl::wait(); }
//从文件的绝对路径中分离文件名并从输出框输出
//从文件的绝对路径中分离文件名并从输出框输出
if (chooser.value() != NULL)
{
strcpy(out_msh_filename,chooser.value());
@ -170,17 +170,17 @@ void cb_mod_file_out_btn(Fl_Widget*, void*){
}
void cb_build_mod_btn(Fl_Widget*, void*){
//首先构建一个GM3D实例
//首先构建一个GM3D实例
GM3D gm3d_instance;
//设置初始化变量
//设置初始化变量
char dimension[1024] = "10/20/990/10/20/990/10/20/490";
char out_mshname[1024] = "Untitled.msh";
char elename[1024] = "Untitled";
//将mesh_para_output拷贝到dimension
//将mesh_para_output拷贝到dimension
strcpy(dimension, mesh_filename);
//std::cout << dimension << std::endl;
//将mod_para_brw中的模型参数列表拷贝到一个临时参数列表
//将mod_para_brw中的模型参数列表拷贝到一个临时参数列表
std::stringstream temp_ss;
modelist temp_list;
modelistArray brw_model_list;
@ -196,21 +196,21 @@ void cb_build_mod_btn(Fl_Widget*, void*){
}
}
//将brw_model_list拷贝到GM3D中
//将brw_model_list拷贝到GM3D中
gm3d_instance.get_model_list(brw_model_list);
//将mod_out_file_input拷贝到out_mshname
//将mod_out_file_input拷贝到out_mshname
strcpy(out_mshname, out_msh_filename);
if (!strcmp(out_mshname,"")){
fl_message("Output file's name can't be empty !");
return;
}
//将mod_ele_input_build拷贝到elename
//将mod_ele_input_build拷贝到elename
strcpy(elename, mod_ele_input_build->value());
if (!strcmp(elename,"")){
fl_message("Element data's name can't be empty !");
return;
}
//构建模型网络
//构建模型网络
if(gm3d_instance.BuildRegularGrid(dimension)){
fl_message("Mesh Parameters Load Error !");
return;
@ -219,12 +219,12 @@ void cb_build_mod_btn(Fl_Widget*, void*){
gm3d_instance.RegisteredOuput(remove_null);
if (gm3d_instance.OutMshFile(out_mshname,elename)) return;
//这里需要一个输出信息窗口
//这里需要一个输出信息窗口
fl_message("Model Construction Completed !");
return;
}
/***********************************以下是正演计算GUI函数******************************/
/***********************************以下是正演计算GUI函数******************************/
char in_msh_filename[1024];
char in_obs_filename[1024];
char out_res_filename[1024];
@ -238,13 +238,13 @@ void cb_mod_file_input(Fl_Input*, void*){
void cb_mod_file_btn(Fl_Widget*, void*){
// Create the file chooser, and show it
// directory, filter, chooser type, title
Fl_File_Chooser chooser(".", "*.msh", Fl_File_Chooser::SINGLE, "Create a model file.");
Fl_File_Chooser chooser(".", "*.msh, *.mst", Fl_File_Chooser::SINGLE, "Create a model file.");
chooser.show();
// Block until user picks something.
while(chooser.shown()) { Fl::wait(); }
//从文件的绝对路径中分离文件名并从输出框输出
//从文件的绝对路径中分离文件名并从输出框输出
if (chooser.value() != NULL)
{
strcpy(in_msh_filename,chooser.value());
@ -279,7 +279,7 @@ void cb_obs_file_btn(Fl_Widget*, void*){
// Block until user picks something.
while(chooser.shown()) { Fl::wait(); }
//从文件的绝对路径中分离文件名并从输出框输出
//从文件的绝对路径中分离文件名并从输出框输出
if (chooser.value() != NULL)
{
strcpy(in_obs_filename,chooser.value());
@ -314,7 +314,7 @@ void cb_res_file_btn(Fl_Widget*, void*){
// Block until user picks something.
while(chooser.shown()) { Fl::wait(); }
//从文件的绝对路径中分离文件名并从输出框输出
//从文件的绝对路径中分离文件名并从输出框输出
if (chooser.value() != NULL)
{
strcpy(out_res_filename,chooser.value());
@ -369,9 +369,9 @@ void cb_mag_data_check(Fl_Check_Button*, void*){
}
void cb_cal_btn(Fl_Button*, void*){
//首先构建一个GM3D实例
//首先构建一个GM3D实例
GM3D gm3d_instance;
//设置初始化变量
//设置初始化变量
char in_mshname[1024] = "Untitled.msh";
char in_obspara[1024] = "Untitled.txt";
char res_outfile[1024] = "Untitled";
@ -386,12 +386,19 @@ void cb_cal_btn(Fl_Button*, void*){
strcpy(ele_name,mod_ele_input->value());
strcpy(noise_para,noise_para_input->value());
strcpy(mag_para,mag_para_input->value());
if(in_mshname[strlen(in_mshname)-1] == 'h'){
if (gm3d_instance.ReadModel(in_mshname,ele_name))
{
fl_message("Model Load Error !");
return;
}
}else if(in_mshname[strlen(in_mshname)-1] == 't'){
if (gm3d_instance.ReadModel_mst(in_mshname))
{
fl_message("Model Load Error !");
return;
}
}
if (gm3d_instance.InitObs(in_obspara))
{
@ -534,6 +541,7 @@ void cb_cal_btn(Fl_Button*, void*){
}
}
fl_message("Forward calculation completed !");
return;
}

30
src/head_func.cpp
View File

@ -1,6 +1,6 @@
#include "head_func.h"
/*************************数据结构函数********************************/
//cpoint减法
/*************************数据结构函数********************************/
//cpoint减法
cpoint operator -(cpoint a, cpoint b){
cpoint m;
m.x=a.x-b.x;
@ -8,25 +8,25 @@ cpoint operator -(cpoint a, cpoint b){
m.z=a.z-b.z;
return m;
}
//cpoint模长
//cpoint模长
double modCpoint(cpoint v){
return sqrt(v.x*v.x+v.y*v.y+v.z*v.z);
}
/*************************全局函数********************************/
//正负分离的atan函数 正数返回atan 负数返回atan+pi
/*************************全局函数********************************/
//正负分离的atan函数 正数返回atan 负数返回atan+pi
double arctg(double v){
double ang;
if(v>=0) ang=atan(v);
else if(v<0) ang=atan(v)+Pi;
return ang;
}
//将string转换为stringstream
//将string转换为stringstream
stringstream str2ss(string s){
stringstream sstr;
sstr.str(""); sstr.clear(); sstr.str(s);
return sstr;
}
//返回一个cpoint的位置id字符串
//返回一个cpoint的位置id字符串
string cpoint_id(cpoint c){
string vert_id, mid_id;
stringstream sstemp;
@ -40,7 +40,7 @@ string cpoint_id(cpoint c){
vert_id = vert_id + " " + mid_id;
return vert_id;
}
//测试打开输入文件 如果成功则返回0并输出信息 否则返回1
//测试打开输入文件 如果成功则返回0并输出信息 否则返回1
int open_infile(ifstream &infile,char* filename){
infile.open(filename);
if (!infile){
@ -49,7 +49,7 @@ int open_infile(ifstream &infile,char* filename){
}
return 0;
}
//测试打开输出文件 如果成功则返回0并输出信息 否则返回1
//测试打开输出文件 如果成功则返回0并输出信息 否则返回1
int open_outfile(ofstream &outfile,char* filename){
outfile.open(filename);
if (!outfile){
@ -59,8 +59,8 @@ int open_outfile(ofstream &outfile,char* filename){
return 0;
}
//规则网络插值 长方形内数据插值 距离平方反比
/*长方体示意图*/
//规则网络插值 长方形内数据插值 距离平方反比
/*长方体示意图*/
// y
// |
// |
@ -69,10 +69,10 @@ int open_outfile(ofstream &outfile,char* filename){
// | |
// | |
// 0------------1--->x
// 左下角坐标x0 y0
// 块体尺寸dx dy
// 插值点坐标x y
// 四个角点值
// 左下角坐标x0 y0
// 块体尺寸dx dy
// 插值点坐标x y
// 四个角点值
double grid_interpolate(double x0,double y0,double dx,double dy,double x,double y,
double d0,double d1,double d2,double d3)
{

98
src/head_func.h
View File

@ -8,7 +8,7 @@
#include "iomanip"
#include "stdio.h"
#include "stdlib.h"
#include "unistd.h"
//#include "unistd.h"
#include "vector"
#include "map"
#include "algorithm"
@ -18,42 +18,42 @@
using namespace std;
//数学常量
#define BDL_MAX 1e+30 ///< 定义变量的最大值
#define BDL_MIN -1e+30 ///< 定义变量的最小值
#define PRECISION 16 ///< 定义小数点后需要使用的位数
#define ZERO 1e-16 ///< 定义零值
//数学常量
#define BDL_MAX 1e+30 ///< 定义变量的最大值
#define BDL_MIN -1e+30 ///< 定义变量的最小值
#define PRECISION 16 ///< 定义小数点后需要使用的位数
#define ZERO 1e-16 ///< 定义零值
//物理常量
#define Pi (4.0*atan(1.0)) ///< 圆周率
#define G0 6.67408e-3 ///< 万有引力常数。注意这里的量级本来应该是e-11，考虑到单位转换，取维度单位为m，密度单位为g/cm^3，乘以G0则重力单位即为mGal
#define T0 5.0e+4 ///< 地磁场平均强度
//宏函数
#define MAX(a,b) (a>b?a:b) ///< 返回a与b的最大值
#define MIN(a,b) (a<b?a:b) ///< 返回a与b的最小值
#define SetToBox(a,b,in) (MAX(a,MIN(b,in))) ///< 返回a与b之间的一个值，若in在a与b之间则直接返回，否则返回较近的边界值
//终端显示控制符
#define BOLDRED "\033[1m\033[31m" ///< 设置后续字符字体为红色加粗
#define BOLDGREEN "\033[1m\033[32m" ///< 设置后续字符字体为绿色加粗
#define BOLDYELLOW "\033[1m\033[33m" ///< 设置后续字符字体为黄色加粗
#define BOLDBLUE "\033[1m\033[34m" ///< 设置后续字符字体为蓝色加粗
#define UNDERLINE "\033[1m\033[4m" ///< 设置后续字符为添加下划线
#define RESET "\033[0m" ///< 重置字符设置
#define MOVEUP(x) printf("\033[%dA", (x)) ///< 将光标向上挪x行
#define MOVEDOWN(x) printf("\033[%dB", (x)) ///< 将光标向下娜x行
#define MOVELEFT(x) printf("\033[%dD", (x)) ///< 将光标向左娜x字符
#define MOVERIGHT(x) printf("\033[%dC", (x)) ///< 将光标向右娜x字符
#define MOVETO(y,x) printf("\033[%d;%dH", (y), (x)) ///< 将光标向右娜动y字符,向上挪动x字符
#define CLEARLINE "\033[K" ///< 清除本行
#define CLEARALL "\033[2J" ///< 清除终端满屏
//数据结构
typedef vector<int> _1iArray; ///< 整形一维向量
typedef vector<double> _1dArray; ///< 双精度浮点一维向量
typedef vector<string> _1sArray; ///< 字符串一维向量
typedef vector<vector<int> > _2iArray; ///< 整形浮点二维向量
typedef vector<vector<double> > _2dArray; ///< 双精度浮点二维向量
typedef map<int,int> _int2intMap; ///< 整型到整形的映射
//模型块体参数
//物理常量
#define Pi (4.0*atan(1.0)) ///< 圆周率
#define G0 6.67408e-3 ///< 万有引力常数。注意这里的量级本来应该是e-11，考虑到单位转换，取维度单位为m，密度单位为g/cm^3，乘以G0则重力单位即为mGal
#define T0 5.0e+4 ///< 地磁场平均强度
//宏函数
#define MAX(a,b) (a>b?a:b) ///< 返回a与b的最大值
#define MIN(a,b) (a<b?a:b) ///< 返回a与b的最小值
#define SetToBox(a,b,in) (MAX(a,MIN(b,in))) ///< 返回a与b之间的一个值，若in在a与b之间则直接返回，否则返回较近的边界值
//终端显示控制符
#define BOLDRED "\033[1m\033[31m" ///< 设置后续字符字体为红色加粗
#define BOLDGREEN "\033[1m\033[32m" ///< 设置后续字符字体为绿色加粗
#define BOLDYELLOW "\033[1m\033[33m" ///< 设置后续字符字体为黄色加粗
#define BOLDBLUE "\033[1m\033[34m" ///< 设置后续字符字体为蓝色加粗
#define UNDERLINE "\033[1m\033[4m" ///< 设置后续字符为添加下划线
#define RESET "\033[0m" ///< 重置字符设置
#define MOVEUP(x) printf("\033[%dA", (x)) ///< 将光标向上挪x行
#define MOVEDOWN(x) printf("\033[%dB", (x)) ///< 将光标向下娜x行
#define MOVELEFT(x) printf("\033[%dD", (x)) ///< 将光标向左娜x字符
#define MOVERIGHT(x) printf("\033[%dC", (x)) ///< 将光标向右娜x字符
#define MOVETO(y,x) printf("\033[%d;%dH", (y), (x)) ///< 将光标向右娜动y字符,向上挪动x字符
#define CLEARLINE "\033[K" ///< 清除本行
#define CLEARALL "\033[2J" ///< 清除终端满屏
//数据结构
typedef vector<int> _1iArray; ///< 整形一维向量
typedef vector<double> _1dArray; ///< 双精度浮点一维向量
typedef vector<string> _1sArray; ///< 字符串一维向量
typedef vector<vector<int> > _2iArray; ///< 整形浮点二维向量
typedef vector<vector<double> > _2dArray; ///< 双精度浮点二维向量
typedef map<int,int> _int2intMap; ///< 整型到整形的映射
//模型块体参数
struct modelist{
char mod_type[1024];
char val_type[1024];
@ -61,33 +61,33 @@ struct modelist{
double mod_value;
};
typedef vector<modelist> modelistArray;
//直角坐标系点
//直角坐标系点
struct cpoint{
int id = -1;
double x = BDL_MAX; double y = BDL_MAX; double z = BDL_MAX;
};
typedef vector<cpoint> cpointArray;
typedef map<string,cpoint> _str2pointMap;
//观测点
//观测点
struct obspoint : public cpoint{
double val = BDL_MAX; double dev = BDL_MAX;
};
typedef vector<obspoint> obspointArray;
//块体
//块体
struct cube{
cpoint cen;
int ids[8] = {-1,-1,-1,-1,-1,-1,-1,-1};
double dx = BDL_MAX; double dy = BDL_MAX; double dz = BDL_MAX;
};
typedef vector<cube> cubeArray;
/*************************数据结构函数********************************/
cpoint operator -(cpoint,cpoint); //矢量减法
double modCpoint(cpoint); //矢量模
/*************************全局函数********************************/
double arctg(double); //正负分离的atan函数 正数返回atan 负数返回atan+pi
stringstream str2ss(string); //将string转换为stringstream
string cpoint_id(cpoint); //返回一个cpoint的位置id
int open_infile(ifstream&,char*); //测试打开输入文件 如果成功则返回0并输出信息 否则返回1
int open_outfile(ofstream&,char*); //测试打开输出文件 如果成功则返回0并输出信息 否则返回1
double grid_interpolate(double,double,double,double,double,double,double,double,double,double); //规则网络插值
/*************************数据结构函数********************************/
cpoint operator -(cpoint,cpoint); //矢量减法
double modCpoint(cpoint); //矢量模
/*************************全局函数********************************/
double arctg(double); //正负分离的atan函数 正数返回atan 负数返回atan+pi
stringstream str2ss(string); //将string转换为stringstream
string cpoint_id(cpoint); //返回一个cpoint的位置id
int open_infile(ifstream&,char*); //测试打开输入文件 如果成功则返回0并输出信息 否则返回1
int open_outfile(ofstream&,char*); //测试打开输出文件 如果成功则返回0并输出信息 否则返回1
double grid_interpolate(double,double,double,double,double,double,double,double,double,double); //规则网络插值
#endif

6
src/init_obs.cpp
View File

@ -8,7 +8,7 @@ int GM3D::InitObs(char* obs_para){
double xmin,xmax,ymin,ymax;
double xs,xe,ys,ye,eleva,dx,dy;
//按格式解析参数 初始化观测位置 用于正演计算
//按格式解析参数 初始化观测位置 用于正演计算
if (7 == sscanf(obs_para,"%lf/%lf/%lf/%lf/%lf/%lf/%lf",&xs,&dx,&xe,&ys,&dy,&ye,&eleva)){
xmin = MIN(xs,xe); xmax = MAX(xs,xe);
ymin = MIN(ys,ye); ymax = MAX(ys,ye);
@ -26,7 +26,7 @@ int GM3D::InitObs(char* obs_para){
x += dx;
}
}
//解析失败 按文件读入 用于反演使用或者正演计算
//解析失败 按文件读入 用于反演使用或者正演计算
else{
ifstream infile;
if (open_infile(infile,obs_para)) return -1;
@ -34,7 +34,7 @@ int GM3D::InitObs(char* obs_para){
while(getline(infile,temp_str)){
if (*(temp_str.begin()) == '#') continue;
else{
//按每行3个数据解析 初始化为用于正演的观测点
//按每行3个数据解析 初始化为用于正演的观测点
if (3 == sscanf(temp_str.c_str(),"%lf %lf %lf",&temp_obs.y,&temp_obs.x,&temp_obs.z)){
temp_obs.z *= -1.0;
temp_obs.id = obs_p_.size();

4
src/out_msh_file.cpp
View File

@ -6,9 +6,9 @@ int GM3D::OutMshFile(char* filename,string data_name){
ofstream outfile;
if (open_outfile(outfile,filename)) return -1;
//好啦 我们这里输出的模型类型应该是块体
//好啦 我们这里输出的模型类型应该是块体
outfile<<"$MeshFormat"<<endl<<"2.2 0 8"<<endl<<"$EndMeshFormat"<<endl<<"$Nodes"<<endl<< out_vert_ids_.size() <<endl;
//输出顶点 最后一个不输出
//输出顶点 最后一个不输出
for (int i = 0; i < out_vert_ids_.size(); i++){
outfile << i << " "
<< setprecision(16) << model_vert_[out_vert_ids_[i]].x << " "

18
src/out_neighbor_file.cpp
View File

@ -3,22 +3,22 @@
int GM3D::OutNeighborFile(char* v_name,char* b_name){
/*
if (strcmp(v_name,"NULL") || strcmp(b_name,"NULL")){
//整理块体间的相邻关系 先初始化顶点相邻数组
//整理块体间的相邻关系 先初始化顶点相邻数组
model_vert_neighbor_.resize(vert_num_);
for (int i = 0; i < vert_num_; i++)
model_vert_neighbor_[i].resize(8,-1);
//遍历所有块体整理顶点相邻关系
//遍历所有块体整理顶点相邻关系
for (int i = 0; i < model_num_; i++){
for (int j = 0; j < 8; j++){
model_vert_neighbor_[model_cube_[i].ids[j]][j] = model_cube_[i].cen.id;
}
}
//遍历所有顶点相邻关系 所有共点的块体都被认定为相邻块体
//遍历所有顶点相邻关系 所有共点的块体都被认定为相邻块体
model_cube_neighbor_.resize(model_num_);
//循环顶点相邻列表 都不为-1则相互添加
//循环顶点相邻列表 都不为-1则相互添加
for (int i = 0; i < vert_num_; i++){
for (int n = 0; n < 8; n++){
if (model_vert_neighbor_[i][n] != -1){
@ -32,14 +32,14 @@ int GM3D::OutNeighborFile(char* v_name,char* b_name){
}
}
vector<int>::iterator pos; //整型向量的迭代器
vector<int>::iterator pos; //整型向量的迭代器
for (int i = 0; i < model_num_; i++){
sort(model_cube_neighbor_[i].begin(),model_cube_neighbor_[i].end()); //对顶点序列由小到大排序
pos = unique(model_cube_neighbor_[i].begin(),model_cube_neighbor_[i].end()); //获取重复序列开始的位置
model_cube_neighbor_[i].erase(pos,model_cube_neighbor_[i].end()); //删除重复点
sort(model_cube_neighbor_[i].begin(),model_cube_neighbor_[i].end()); //对顶点序列由小到大排序
pos = unique(model_cube_neighbor_[i].begin(),model_cube_neighbor_[i].end()); //获取重复序列开始的位置
model_cube_neighbor_[i].erase(pos,model_cube_neighbor_[i].end()); //删除重复点
}
//清理数组
//清理数组
for (int i = 0; i < vert_num_; i++){
model_vert_neighbor_[i].clear();
vector <int>().swap(model_vert_neighbor_[i]);

2
src/progress_bar.cpp
View File

@ -47,7 +47,7 @@ int ProgressBar::GetConsoleWidth(){
width = csbi.srWindow.Right - csbi.srWindow.Left;
#else
struct winsize win;
//注意！当我们使用pipe here-doc等通道获取程序参数时无法正确的获取窗口大小 此时我们将使用预定值
//注意！当我们使用pipe here-doc等通道获取程序参数时无法正确的获取窗口大小 此时我们将使用预定值
if (ioctl(0, TIOCGWINSZ, &win) != -1)
width = win.ws_col;
else width = 100;

2
src/progress_bar.h
View File

@ -1,6 +1,6 @@
#ifndef _PROGRESS_BAR_
#define _PROGRESS_BAR_
#include <sys/ioctl.h>
//#include </*sys/ioctl.h*/>
#include <iostream>
#include <iomanip>
#include <cstring>

44
src/read_model.cpp
View File

@ -10,14 +10,14 @@ int GM3D::ReadModel(char* filename,char* input_forward_model_name){
stringstream temp_ss;
ifstream mshin;
if (open_infile(mshin,filename)) return -1; //检查并打开模型文件
if (open_infile(mshin,filename)) return -1; //检查并打开模型文件
while(getline(mshin,temp_str)){
//读入模型空间顶点集 msh文件版本为2.2
//读入模型空间顶点集 msh文件版本为2.2
if (temp_str == "$Nodes"){
getline(mshin,temp_str);
temp_ss = str2ss(temp_str);
temp_ss >> vert_num_; //第一个数为顶点的个数
model_vert_.resize(vert_num_); //开辟空间
temp_ss >> vert_num_; //第一个数为顶点的个数
model_vert_.resize(vert_num_); //开辟空间
for (int i = 0; i < vert_num_; i++){
getline(mshin,temp_str);
temp_ss = str2ss(temp_str);
@ -25,19 +25,19 @@ int GM3D::ReadModel(char* filename,char* input_forward_model_name){
model_vert_[i] = temp_vert;
}
}
//读入模型空间单元体
//读入模型空间单元体
else if (temp_str == "$Elements"){
getline(mshin,temp_str);
temp_ss = str2ss(temp_str);
temp_ss >> model_num_; //第一个数为总元素的个数 包含了所有类型的元素 比如三角形 四边形 块体等
temp_ss >> model_num_; //第一个数为总元素的个数 包含了所有类型的元素 比如三角形 四边形 块体等
model_cube_.resize(model_num_);
for (int i = 0; i < model_num_; i++){
getline(mshin,temp_str);
temp_ss = str2ss(temp_str);
temp_ss >> temp_cu.cen.id >> ele_type;
//只读入块体
//只读入块体
if (ele_type == 5){
temp_ss >> attri_num; //跳过模型单元的几何组与物理组等信息 以后可能会有用
temp_ss >> attri_num; //跳过模型单元的几何组与物理组等信息 以后可能会有用
for (int a = 0; a < attri_num; a++)
temp_ss >> temp_attri;
for (int a = 0; a < 8; a++)
@ -46,44 +46,44 @@ int GM3D::ReadModel(char* filename,char* input_forward_model_name){
}
}
}
else continue; //不能识别的单元都被忽略了
else continue; //不能识别的单元都被忽略了
}
mshin.close();
//第二次读入模型文件 初始化模型单元属性
if (open_infile(mshin,filename)) return -1; //检查并打开模型文件
//第二次读入模型文件 初始化模型单元属性
if (open_infile(mshin,filename)) return -1; //检查并打开模型文件
while(getline(mshin,temp_str)){
//读入模型单元属性 注意因为msh文件中$ElementData并未注明所属元素类型
//所以可能会将其他元素类型的属性值也读入 但因为其在pyIdMap中并未注册 所以属性值会全为0 在后续使用时我们需要通过名称辨别
//读入模型单元属性 注意因为msh文件中$ElementData并未注明所属元素类型
//所以可能会将其他元素类型的属性值也读入 但因为其在pyIdMap中并未注册 所以属性值会全为0 在后续使用时我们需要通过名称辨别
if (temp_str == "$ElementData"){
temp_model.resize(model_num_,0.0); //初始化temp_model 为读入模型单元属性做准备
for (int i = 0; i < 2; i++) //先读入元素块的名称 添加到数组
temp_model.resize(model_num_,0.0); //初始化temp_model 为读入模型单元属性做准备
for (int i = 0; i < 2; i++) //先读入元素块的名称 添加到数组
getline(mshin,temp_str);
input_model_names_.push_back(temp_str);
for (int i = 0; i < 6; i++) //跳过元素属性前面的值 最后一次为当前元素块的个数
for (int i = 0; i < 6; i++) //跳过元素属性前面的值 最后一次为当前元素块的个数
getline(mshin,temp_str);
temp_ss = str2ss(temp_str);
temp_ss >> temp_int;
for (int i = 0; i < temp_int; i++){
getline(mshin,temp_str);
temp_ss = str2ss(temp_str);
temp_ss >> temp_id >> temp_val; //读入单元体索引与属性值
temp_ss >> temp_id >> temp_val; //读入单元体索引与属性值
temp_model[temp_id] = temp_val;
}
input_models_.push_back(temp_model);
temp_model.clear();
}
else continue; //不能识别的单元都被忽略了
else continue; //不能识别的单元都被忽略了
}
mshin.close();
//清理映射
//清理映射
temp_model.clear();
vector <double>().swap(temp_model);
//初始化数组
//初始化数组
forward_model_.resize(model_num_,0.0);
//匹配数据名称
//匹配数据名称
char forward_model_name[1024] = "\"";
strcat(forward_model_name,input_forward_model_name);
strcat(forward_model_name,"\"");
@ -94,7 +94,7 @@ int GM3D::ReadModel(char* filename,char* input_forward_model_name){
}
}
//计算块体的中心位置和尺寸
//计算块体的中心位置和尺寸
cpoint corner[8];
for (int i = 0; i < model_num_; i++){
for (int j = 0; j < 8; j++){

51
src/read_model_mst.cpp Normal file
View File

@ -0,0 +1,51 @@
#include "gm3d.h"
int GM3D::ReadModel_mst(char* filename){
int temp_int,ele_type,attri_num,temp_attri,temp_id;
double temp_val;
_1dArray temp_model;
cpoint temp_vert;
cube temp_cu;
string temp_str;
stringstream temp_ss;
ifstream mstin;
if (open_infile(mstin,filename)) return -1; //检查并打开模型文件
double ynum, xnum, znum;
double ystart, xstart, zstart;
double ylength, xlength, zlength;
char temp;
mstin>>ynum>>xnum>>znum;
mstin>>ystart>>xstart>>zstart;
mstin>>ynum>>temp>>ylength;
mstin>>xnum>>temp>>xlength;
mstin>>znum>>temp>>zlength;
model_num_ = ynum*xnum*znum;
model_cube_.resize(model_num_);
//初始化数组
forward_model_.resize(model_num_,0.0);
double f_mdl_tmp;
int i, j, k;
int tmp_num;
for(i = 0; i < xnum; i++){
for(j = 0; j < ynum; j++){
for(k = 0; k < znum; k++){
tmp_num = i*ynum*znum + j*znum + k;
mstin>>f_mdl_tmp;
forward_model_[tmp_num] = f_mdl_tmp;
model_cube_[tmp_num].cen.x = xstart + xlength*i + xlength*0.5;
model_cube_[tmp_num].cen.y = ystart + ylength*j + ylength*0.5;
model_cube_[tmp_num].cen.z = zstart + zlength*k + zlength*0.5;
model_cube_[tmp_num].dx = xlength;
model_cube_[tmp_num].dy = ylength;
model_cube_[tmp_num].dz = zlength;
}
}
}
mstin.close();
return 0;
}

26
src/registered_output.cpp
View File

@ -2,10 +2,10 @@
int GM3D::RegisteredOuput(bool remove_empty_element){
int count;
//统计输出模型单元块体和顶点 以及输出的块体数据列表
//统计输出模型单元块体和顶点 以及输出的块体数据列表
if (remove_empty_element){
count = 0;
//遍历所有块体数据 注册有值的块体
//遍历所有块体数据 注册有值的块体
for (int i = 0; i < model_num_; i++){
if (model_block_val_[i] != BDL_MAX){
out_ele_data_ids_.push_back(i);
@ -15,38 +15,38 @@ int GM3D::RegisteredOuput(bool remove_empty_element){
}
}
//遍历所有注册的块体 添加顶点
//遍历所有注册的块体 添加顶点
for (int i = 0; i < out_ele_ids_.size(); i++){
for (int j = 0; j < 8; j++){
out_vert_ids_.push_back(model_cube_[out_ele_ids_[i]].ids[j]);
}
}
//去除输出顶点中的重复部分
vector<int>::iterator pos; //整型向量的迭代器
sort(out_vert_ids_.begin(),out_vert_ids_.end()); //对顶点序列由小到大排序
pos = unique(out_vert_ids_.begin(),out_vert_ids_.end()); //获取重复序列开始的位置
out_vert_ids_.erase(pos,out_vert_ids_.end()); //删除重复点
//去除输出顶点中的重复部分
vector<int>::iterator pos; //整型向量的迭代器
sort(out_vert_ids_.begin(),out_vert_ids_.end()); //对顶点序列由小到大排序
pos = unique(out_vert_ids_.begin(),out_vert_ids_.end()); //获取重复序列开始的位置
out_vert_ids_.erase(pos,out_vert_ids_.end()); //删除重复点
//将需要输出的模型顶点序号与它们的排序做一个对应 保证在输出文件中顶点索引号始终是从0开始的连续的序列
//将需要输出的模型顶点序号与它们的排序做一个对应 保证在输出文件中顶点索引号始终是从0开始的连续的序列
for (int i = 0; i < out_vert_ids_.size(); i++){
vert_out_map_[out_vert_ids_[i]] = i;
}
}
else{
//输出的模型块体为所有
//输出的模型块体为所有
out_ele_ids_.resize(model_num_);
for (int i = 0; i < model_num_; i++){
out_ele_ids_[i] = i;
}
//输出的模型顶点为所有
//输出的模型顶点为所有
out_vert_ids_.resize(vert_num_);
for (int i = 0; i < vert_num_; i++){
vert_out_map_[i] = out_vert_ids_[i] = i;
}
//注册所有有值的块体数据
//注册所有有值的块体数据
count = 0;
for (int i = 0; i < model_num_; i++){
if (model_block_val_[i] != BDL_MAX){
@ -56,6 +56,6 @@ int GM3D::RegisteredOuput(bool remove_empty_element){
}
}
}
//输出所有模型单元块体和顶点 只统计输出的块体数据列表
//输出所有模型单元块体和顶点 只统计输出的块体数据列表
return 0;
}