zhangyiss/separation-of-potential-data
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

tmp update

Browse Source
This commit is contained in:
张壹
2022-02-08 17:29:12 +08:00
parent a95cc8a31c
commit 7254ab3fe8
148 changed files with 624322 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

17
src/CMakeLists.txt Normal file
View File

@@ -0,0 +1,17 @@
aux_source_directory(. SRC_DIR)
set(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/bin)
add_executable(sgm_esl ${SRC_DIR})
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --std=c++11 -O3")
# 添加openmp的编译命令 设置编译选项
find_package(OpenMP REQUIRED)
if (OpenMP_CXX_FOUND)
message(STATUS "OpenMP Found.")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}")
set(CMAKE_SHARED_LINKER_FLAGS "${CMAKE_SHARED_LINKER_FLAGS} ${OpenMP_SHARED_LINKER_FLAGS}")
target_link_libraries(sgm_esl PUBLIC OpenMP::OpenMP_CXX)
endif()

137
src/cal_deltaT_kernel.cpp Normal file
View File

@@ -0,0 +1,137 @@
#include "sgm_esl.h"
#include "progressBar.h"
void SGM_ESL::CalDeltaTKernel(){
int i,j;
//double k1,k2,k3,k4,k5,k6;
double x1,x2,y1,y2,z1,z2;
double w111,w112,w121,w122,w211,w212,w221,w222;
double v111,v112,v121,v122,v211,v212,v221,v222;
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
double *k1, *k2, *k3, *k4, *k5, *k6;
k1 = new double [mod_blocks_num_];
k2 = new double [mod_blocks_num_];
k3 = new double [mod_blocks_num_];
k4 = new double [mod_blocks_num_];
k5 = new double [mod_blocks_num_];
k6 = new double [mod_blocks_num_];
double I, D;
I0_ = I0_ * Pi / 180.0;
D0_ = D0_ * Pi / 180.0;
// 如果区域属性被激活 则使用区域磁化参数 否则使用局部磁化参数
for (int i = 0; i < mod_blocks_num_; i++)
{
if (regional_model_list_[i] == 1.0)
{
I = reg_I_ * Pi / 180.0;
D = reg_D_ * Pi / 180.0;
}
else
{
I = loc_I_ * Pi / 180.0;
D = loc_D_ * Pi / 180.0;
}
k1[i]=cos(I0_)*sin(D0_)*sin(I)+sin(I0_)*cos(I)*sin(D);
k2[i]=cos(I0_)*cos(D0_)*sin(I)+sin(I0_)*cos(I)*cos(D);
k3[i]=cos(I0_)*cos(D0_)*cos(I)*sin(D)+cos(I0_)*sin(D0_)*cos(I)*cos(D);
k4[i]=cos(I0_)*cos(D0_)*cos(I)*cos(D);
k5[i]=cos(I0_)*sin(D0_)*cos(I)*sin(D);
k6[i]=-sin(I0_)*sin(I);
}
if (GM_kernel_.empty())
{
GM_kernel_.resize(obs_points_num_);
for (i = 0; i < obs_points_num_; i++){
GM_kernel_[i].resize(mod_blocks_num_,0.0);
}
}
ProgressBar *bar = new ProgressBar(obs_points_num_,"Calculating GMKernel");
for (i = 0; i < obs_points_num_; i++){
bar->Progressed(i);
#pragma omp parallel for private(j,R222,R122,R212,R112,R221,R121,R211,R111,\
G222,G122,G212,G112,G221,G121,G211,G111,\
w111,w112,w121,w122,w211,w212,w221,w222,\
v111,v112,v121,v122,v211,v212,v221,v222,\
x1,x2,y1,y2,z1,z2) shared(i) schedule(guided)
for (j = 0; j < mod_blocks_num_; j++){
x1 = mod_blocks_[j].xmin; x2 = mod_blocks_[j].xmax;
y1 = mod_blocks_[j].ymin; y2 = mod_blocks_[j].ymax;
z1 = mod_blocks_[j].zmin; z2 = mod_blocks_[j].zmax;
R222=sqrt((x2-obs_points_[i].x)*(x2-obs_points_[i].x)+(y2-obs_points_[i].y)*(y2-obs_points_[i].y)+(z2-obs_points_[i].z)*(z2-obs_points_[i].z));
R122=sqrt((x1-obs_points_[i].x)*(x1-obs_points_[i].x)+(y2-obs_points_[i].y)*(y2-obs_points_[i].y)+(z2-obs_points_[i].z)*(z2-obs_points_[i].z));
R212=sqrt((x2-obs_points_[i].x)*(x2-obs_points_[i].x)+(y1-obs_points_[i].y)*(y1-obs_points_[i].y)+(z2-obs_points_[i].z)*(z2-obs_points_[i].z));
R112=sqrt((x1-obs_points_[i].x)*(x1-obs_points_[i].x)+(y1-obs_points_[i].y)*(y1-obs_points_[i].y)+(z2-obs_points_[i].z)*(z2-obs_points_[i].z));
R221=sqrt((x2-obs_points_[i].x)*(x2-obs_points_[i].x)+(y2-obs_points_[i].y)*(y2-obs_points_[i].y)+(z1-obs_points_[i].z)*(z1-obs_points_[i].z));
R121=sqrt((x1-obs_points_[i].x)*(x1-obs_points_[i].x)+(y2-obs_points_[i].y)*(y2-obs_points_[i].y)+(z1-obs_points_[i].z)*(z1-obs_points_[i].z));
R211=sqrt((x2-obs_points_[i].x)*(x2-obs_points_[i].x)+(y1-obs_points_[i].y)*(y1-obs_points_[i].y)+(z1-obs_points_[i].z)*(z1-obs_points_[i].z));
R111=sqrt((x1-obs_points_[i].x)*(x1-obs_points_[i].x)+(y1-obs_points_[i].y)*(y1-obs_points_[i].y)+(z1-obs_points_[i].z)*(z1-obs_points_[i].z));
w222=((x2-obs_points_[i].x)*(x2-obs_points_[i].x))+R222*(z2-obs_points_[i].z)+((z2-obs_points_[i].z)*(z2-obs_points_[i].z));
w122=((x1-obs_points_[i].x)*(x1-obs_points_[i].x))+R122*(z2-obs_points_[i].z)+((z2-obs_points_[i].z)*(z2-obs_points_[i].z));
w212=((x2-obs_points_[i].x)*(x2-obs_points_[i].x))+R212*(z2-obs_points_[i].z)+((z2-obs_points_[i].z)*(z2-obs_points_[i].z));
w112=((x1-obs_points_[i].x)*(x1-obs_points_[i].x))+R112*(z2-obs_points_[i].z)+((z2-obs_points_[i].z)*(z2-obs_points_[i].z));
w221=((x2-obs_points_[i].x)*(x2-obs_points_[i].x))+R221*(z1-obs_points_[i].z)+((z1-obs_points_[i].z)*(z1-obs_points_[i].z));
w121=((x1-obs_points_[i].x)*(x1-obs_points_[i].x))+R121*(z1-obs_points_[i].z)+((z1-obs_points_[i].z)*(z1-obs_points_[i].z));
w211=((x2-obs_points_[i].x)*(x2-obs_points_[i].x))+R211*(z1-obs_points_[i].z)+((z1-obs_points_[i].z)*(z1-obs_points_[i].z));
w111=((x1-obs_points_[i].x)*(x1-obs_points_[i].x))+R111*(z1-obs_points_[i].z)+((z1-obs_points_[i].z)*(z1-obs_points_[i].z));
v222=((y2-obs_points_[i].y)*(y2-obs_points_[i].y))+R222*(z2-obs_points_[i].z)+((z2-obs_points_[i].z)*(z2-obs_points_[i].z));
v122=((y2-obs_points_[i].y)*(y2-obs_points_[i].y))+R122*(z2-obs_points_[i].z)+((z2-obs_points_[i].z)*(z2-obs_points_[i].z));
v212=((y1-obs_points_[i].y)*(y1-obs_points_[i].y))+R212*(z2-obs_points_[i].z)+((z2-obs_points_[i].z)*(z2-obs_points_[i].z));
v112=((y1-obs_points_[i].y)*(y1-obs_points_[i].y))+R112*(z2-obs_points_[i].z)+((z2-obs_points_[i].z)*(z2-obs_points_[i].z));
v221=((y2-obs_points_[i].y)*(y2-obs_points_[i].y))+R221*(z1-obs_points_[i].z)+((z1-obs_points_[i].z)*(z1-obs_points_[i].z));
v121=((y2-obs_points_[i].y)*(y2-obs_points_[i].y))+R121*(z1-obs_points_[i].z)+((z1-obs_points_[i].z)*(z1-obs_points_[i].z));
v211=((y1-obs_points_[i].y)*(y1-obs_points_[i].y))+R211*(z1-obs_points_[i].z)+((z1-obs_points_[i].z)*(z1-obs_points_[i].z));
v111=((y1-obs_points_[i].y)*(y1-obs_points_[i].y))+R111*(z1-obs_points_[i].z)+((z1-obs_points_[i].z)*(z1-obs_points_[i].z));
G222=k1[j]*log(R222+x2-obs_points_[i].x)+k2[j]*log(R222+y2-obs_points_[i].y)+k3[j]*log(R222+z2-obs_points_[i].z)
+k4[j]*atan(((x2-obs_points_[i].x)*(y2-obs_points_[i].y))/w222)+k5[j]*atan(((x2-obs_points_[i].x)*(y2-obs_points_[i].y))/v222)
+k6[j]*atan(((x2-obs_points_[i].x)*(y2-obs_points_[i].y))/(R222*(z2-obs_points_[i].z)));
G122=k1[j]*log(R122+x1-obs_points_[i].x)+k2[j]*log(R122+y2-obs_points_[i].y)+k3[j]*log(R122+z2-obs_points_[i].z)
+k4[j]*atan(((x1-obs_points_[i].x)*(y2-obs_points_[i].y))/w122)+k5[j]*atan(((x1-obs_points_[i].x)*(y2-obs_points_[i].y))/v122)
+k6[j]*atan(((x1-obs_points_[i].x)*(y2-obs_points_[i].y))/(R122*(z2-obs_points_[i].z)));
G212=k1[j]*log(R212+x2-obs_points_[i].x)+k2[j]*log(R212+y1-obs_points_[i].y)+k3[j]*log(R212+z2-obs_points_[i].z)
+k4[j]*atan(((x2-obs_points_[i].x)*(y1-obs_points_[i].y))/w212)+k5[j]*atan(((x2-obs_points_[i].x)*(y1-obs_points_[i].y))/v212)
+k6[j]*atan(((x2-obs_points_[i].x)*(y1-obs_points_[i].y))/(R212*(z2-obs_points_[i].z)));
G112=k1[j]*log(R112+x1-obs_points_[i].x)+k2[j]*log(R112+y1-obs_points_[i].y)+k3[j]*log(R112+z2-obs_points_[i].z)
+k4[j]*atan(((x1-obs_points_[i].x)*(y1-obs_points_[i].y))/w112)+k5[j]*atan(((x1-obs_points_[i].x)*(y1-obs_points_[i].y))/v112)
+k6[j]*atan(((x1-obs_points_[i].x)*(y1-obs_points_[i].y))/(R112*(z2-obs_points_[i].z)));
G221=k1[j]*log(R221+x2-obs_points_[i].x)+k2[j]*log(R221+y2-obs_points_[i].y)+k3[j]*log(R221+z1-obs_points_[i].z)
+k4[j]*atan(((x2-obs_points_[i].x)*(y2-obs_points_[i].y))/w221)+k5[j]*atan(((x2-obs_points_[i].x)*(y2-obs_points_[i].y))/v221)
+k6[j]*atan(((x2-obs_points_[i].x)*(y2-obs_points_[i].y))/(R221*(z1-obs_points_[i].z)));
G121=k1[j]*log(R121+x1-obs_points_[i].x)+k2[j]*log(R121+y2-obs_points_[i].y)+k3[j]*log(R121+z1-obs_points_[i].z)
+k4[j]*atan(((x1-obs_points_[i].x)*(y2-obs_points_[i].y))/w121)+k5[j]*atan(((x1-obs_points_[i].x)*(y2-obs_points_[i].y))/v121)
+k6[j]*atan(((x1-obs_points_[i].x)*(y2-obs_points_[i].y))/(R121*(z1-obs_points_[i].z)));
G211=k1[j]*log(R211+x2-obs_points_[i].x)+k2[j]*log(R211+y1-obs_points_[i].y)+k3[j]*log(R211+z1-obs_points_[i].z)
+k4[j]*atan(((x2-obs_points_[i].x)*(y1-obs_points_[i].y))/w211)+k5[j]*atan(((x2-obs_points_[i].x)*(y1-obs_points_[i].y))/v211)
+k6[j]*atan(((x2-obs_points_[i].x)*(y1-obs_points_[i].y))/(R211*(z1-obs_points_[i].z)));
G111=k1[j]*log(R111+x1-obs_points_[i].x)+k2[j]*log(R111+y1-obs_points_[i].y)+k3[j]*log(R111+z1-obs_points_[i].z)
+k4[j]*atan(((x1-obs_points_[i].x)*(y1-obs_points_[i].y))/w111)+k5[j]*atan(((x1-obs_points_[i].x)*(y1-obs_points_[i].y))/v111)
+k6[j]*atan(((x1-obs_points_[i].x)*(y1-obs_points_[i].y))/(R111*(z1-obs_points_[i].z)));
GM_kernel_[i][j] = T0/(4*Pi)*(G222-G122-G212+G112-G221+G121+G211-G111);
}
}
delete[] k1;
delete[] k2;
delete[] k3;
delete[] k4;
delete[] k5;
delete[] k6;
return;
}

46
src/cal_g_kernel.cpp Normal file
View File

@@ -0,0 +1,46 @@
#include "sgm_esl.h"
#include "progressBar.h"
void SGM_ESL::CalGKernel(){
int i, j;
double x1,x2,y1,y2,z1,z2;
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
GM_kernel_.resize(obs_points_num_);
for (i = 0; i < obs_points_num_; i++){
GM_kernel_[i].resize(mod_blocks_num_,0.0);
}
ProgressBar *bar = new ProgressBar(obs_points_num_,"Calculating GMKernel");
for (i = 0; i < obs_points_num_; i++){
bar->Progressed(i);
#pragma omp parallel for private(j,x1,x2,y1,y2,z1,z2,R222,R122,R212,R112,R221,R121,R211,R111,G222,G122,G212,G112,G221,G121,G211,G111) shared(i) schedule(guided)
for (j = 0; j < mod_blocks_num_; j++){
x1 = mod_blocks_[j].xmin; x2 = mod_blocks_[j].xmax;
y1 = mod_blocks_[j].ymin; y2 = mod_blocks_[j].ymax;
z1 = mod_blocks_[j].zmin; z2 = mod_blocks_[j].zmax;
R222=sqrt((x2- obs_points_[i].x)*(x2- obs_points_[i].x)+(y2- obs_points_[i].y)*(y2- obs_points_[i].y)+(z2- obs_points_[i].z)*(z2- obs_points_[i].z));
R122=sqrt((x1- obs_points_[i].x)*(x1- obs_points_[i].x)+(y2- obs_points_[i].y)*(y2- obs_points_[i].y)+(z2- obs_points_[i].z)*(z2- obs_points_[i].z));
R212=sqrt((x2- obs_points_[i].x)*(x2- obs_points_[i].x)+(y1- obs_points_[i].y)*(y1- obs_points_[i].y)+(z2- obs_points_[i].z)*(z2- obs_points_[i].z));
R112=sqrt((x1- obs_points_[i].x)*(x1- obs_points_[i].x)+(y1- obs_points_[i].y)*(y1- obs_points_[i].y)+(z2- obs_points_[i].z)*(z2- obs_points_[i].z));
R221=sqrt((x2- obs_points_[i].x)*(x2- obs_points_[i].x)+(y2- obs_points_[i].y)*(y2- obs_points_[i].y)+(z1- obs_points_[i].z)*(z1- obs_points_[i].z));
R121=sqrt((x1- obs_points_[i].x)*(x1- obs_points_[i].x)+(y2- obs_points_[i].y)*(y2- obs_points_[i].y)+(z1- obs_points_[i].z)*(z1- obs_points_[i].z));
R211=sqrt((x2- obs_points_[i].x)*(x2- obs_points_[i].x)+(y1- obs_points_[i].y)*(y1- obs_points_[i].y)+(z1- obs_points_[i].z)*(z1- obs_points_[i].z));
R111=sqrt((x1- obs_points_[i].x)*(x1- obs_points_[i].x)+(y1- obs_points_[i].y)*(y1- obs_points_[i].y)+(z1- obs_points_[i].z)*(z1- obs_points_[i].z));
G222=(x2- obs_points_[i].x)*log((y2- obs_points_[i].y)+R222)+(y2- obs_points_[i].y)*log((x2- obs_points_[i].x)+R222)+(z2- obs_points_[i].z)*arctg((z2- obs_points_[i].z)*R222/(x2- obs_points_[i].x)/(y2- obs_points_[i].y));
G122=(x1- obs_points_[i].x)*log((y2- obs_points_[i].y)+R122)+(y2- obs_points_[i].y)*log((x1- obs_points_[i].x)+R122)+(z2- obs_points_[i].z)*arctg((z2- obs_points_[i].z)*R122/(x1- obs_points_[i].x)/(y2- obs_points_[i].y));
G212=(x2- obs_points_[i].x)*log((y1- obs_points_[i].y)+R212)+(y1- obs_points_[i].y)*log((x2- obs_points_[i].x)+R212)+(z2- obs_points_[i].z)*arctg((z2- obs_points_[i].z)*R212/(x2- obs_points_[i].x)/(y1- obs_points_[i].y));
G112=(x1- obs_points_[i].x)*log((y1- obs_points_[i].y)+R112)+(y1- obs_points_[i].y)*log((x1- obs_points_[i].x)+R112)+(z2- obs_points_[i].z)*arctg((z2- obs_points_[i].z)*R112/(x1- obs_points_[i].x)/(y1- obs_points_[i].y));
G221=(x2- obs_points_[i].x)*log((y2- obs_points_[i].y)+R221)+(y2- obs_points_[i].y)*log((x2- obs_points_[i].x)+R221)+(z1- obs_points_[i].z)*arctg((z1- obs_points_[i].z)*R221/(x2- obs_points_[i].x)/(y2- obs_points_[i].y));
G121=(x1- obs_points_[i].x)*log((y2- obs_points_[i].y)+R121)+(y2- obs_points_[i].y)*log((x1- obs_points_[i].x)+R121)+(z1- obs_points_[i].z)*arctg((z1- obs_points_[i].z)*R121/(x1- obs_points_[i].x)/(y2- obs_points_[i].y));
G211=(x2- obs_points_[i].x)*log((y1- obs_points_[i].y)+R211)+(y1- obs_points_[i].y)*log((x2- obs_points_[i].x)+R211)+(z1- obs_points_[i].z)*arctg((z1- obs_points_[i].z)*R211/(x2- obs_points_[i].x)/(y1- obs_points_[i].y));
G111=(x1- obs_points_[i].x)*log((y1- obs_points_[i].y)+R111)+(y1- obs_points_[i].y)*log((x1- obs_points_[i].x)+R111)+(z1- obs_points_[i].z)*arctg((z1- obs_points_[i].z)*R111/(x1- obs_points_[i].x)/(y1- obs_points_[i].y));
GM_kernel_[i][j] = -1.0*G0*(G222-G122-G212+G112-G221+G121+G211-G111);
}
}
return;
}

17
src/cal_part_b.cpp Normal file
View File

@@ -0,0 +1,17 @@
#include "sgm_esl.h"
void SGM_ESL::CalPartB(_1dArray obs_array){
int i, j;
//PartB = G^T * W_d^TW_d * d_obs + mu * W_m^TW_m * M_ref
//W_m^TW_m = WzTWz*WvTWv*wm
#pragma omp parallel for private(j) schedule(guided)
for (j = 0; j < mod_blocks_num_; j++)
{
PartB[j] = 0.0;
for (i = 0; i < obs_points_num_; i++)
{
PartB[j] += GM_kernel_[i][j]*WdTWd[i]*obs_array[i];
}
}
return;
}

27
src/cal_partial_model.cpp Normal file
View File

@@ -0,0 +1,27 @@
#include "sgm_esl.h"
void SGM_ESL::CalPartialModel(_1dArray& out_obs_array)
{
int i,j;
#pragma omp parallel for private(i,j) schedule(guided)
for (i = 0; i < obs_points_num_; i++){
out_obs_array[i] = 0.0;
for (j = 0; j < mod_blocks_num_; j++){
out_obs_array[i] += GM_kernel_[i][j] * invert_model_[j];
}
}
return;
}
void SGM_ESL::CalPartialModel(_1dArray& out_obs_array,_1dArray weight_list)
{
int i,j;
#pragma omp parallel for private(i,j) schedule(guided)
for (i = 0; i < obs_points_num_; i++){
out_obs_array[i] = 0.0;
for (j = 0; j < mod_blocks_num_; j++){
out_obs_array[i] += GM_kernel_[i][j] * invert_model_[j] * weight_list[j];
}
}
return;
}

15
src/cal_wd.cpp Normal file
View File

@@ -0,0 +1,15 @@
#include "sgm_esl.h"
#include "progressBar.h"
void SGM_ESL::CalWd(double attun_factor){
//reserve vector
if (!WdTWd.empty()) WdTWd.clear();
WdTWd.resize(obs_points_num_,0.0);
ProgressBar *bar = new ProgressBar(obs_points_num_,"Calculating Wd");
for (int i = 0; i < obs_points_num_; i++)
{
bar->Progressed(i);
WdTWd[i] = attun_factor/(obs_points_[i].dev*obs_points_[i].dev*obs_points_num_);
}
return;
}

12
src/cal_wp.cpp Normal file
View File

@@ -0,0 +1,12 @@
#include "sgm_esl.h"
#include "progressBar.h"
void SGM_ESL::CalWp(double wp_factor){
Wp.resize(mod_blocks_num_,0.0);
ProgressBar *bar = new ProgressBar(mod_blocks_num_,"Calculating Wp");
for (int i = 0; i < mod_blocks_num_; i++){
bar->Progressed(i);
Wp[i] = wp_factor/(WzTWz[i]*WvTWv[i]);
}
return;
}

38
src/cal_wv_wz.cpp Normal file
View File

@@ -0,0 +1,38 @@
#include "sgm_esl.h"
#include "progressBar.h"
void SGM_ESL::CalWvWz(){
double weight_angle;
double max_vol = BDL_MIN, max_depth = BDL_MIN;
cpoint blo_center,obs_center;
WvTWv.resize(mod_blocks_num_,0.0); WzTWz.resize(mod_blocks_num_,0.0);
PartB.resize(mod_blocks_num_,0.0);
//计算体积加权矩阵
ProgressBar *bar = new ProgressBar(mod_blocks_num_,"Calculating Wv");
for (int i = 0; i < mod_blocks_num_; i++){
bar->Progressed(i);
WvTWv[i] = (mod_blocks_[i].xmax - mod_blocks_[i].xmin) * (mod_blocks_[i].ymax - mod_blocks_[i].ymin) * (mod_blocks_[i].zmax - mod_blocks_[i].zmin);
max_vol = MAX(max_vol,WvTWv[i]);
}
//计算深度加权矩阵
ProgressBar *bar2 = new ProgressBar(mod_blocks_num_,"Calculating Wz");
for (int i = 0; i < mod_blocks_num_; i++){
bar2->Progressed(i);
blo_center.x = 0.5*(mod_blocks_[i].xmin + mod_blocks_[i].xmax);
blo_center.y = 0.5*(mod_blocks_[i].ymin + mod_blocks_[i].ymax);
blo_center.z = 0.5*(mod_blocks_[i].zmin + mod_blocks_[i].zmax);
for (int j = 0; j < obs_points_num_; j++){
weight_angle = (blo_center.z + obs_points_[j].z)/modCpoint(blo_center - obs_points_[j]);
WzTWz[i] += weight_angle * pow(1.0/pow(modCpoint(obs_points_[j]-blo_center)+ZERO,3),2);
}
max_depth = MAX(max_depth,WzTWz[i]);
}
for (int i = 0; i < mod_blocks_.size(); i++){
WvTWv[i] = pow(WvTWv[i]/max_vol,alpha_);
WzTWz[i] = pow(WzTWz[i]/max_depth,beta_);
}
return;
}

77
src/head_func.cpp Normal file
View File

@@ -0,0 +1,77 @@
#include "head_func.h"
/*************************矢量函数********************************/
//cpoint减法
cpoint operator -(cpoint a, cpoint b){
cpoint m;
m.x=a.x-b.x;
m.y=a.y-b.y;
m.z=a.z-b.z;
return m;
}
//cpoint乘标量 注意标量在前
cpoint operator *(double sign,cpoint b){
cpoint m;
m.x=sign*b.x;
m.y=sign*b.y;
m.z=sign*b.z;
return m;
}
//cpoint模长
double modCpoint(cpoint v){
return sqrt(v.x*v.x+v.y*v.y+v.z*v.z);
}
//cpoint两点向量夹角
double angCpoint(cpoint a,cpoint b)
{
double arg = SetToBox(-1.0,1.0,(a.x*b.x+a.y*b.y+a.z*b.z)/(sqrt(a.x*a.x+a.y*a.y+a.z*a.z)*sqrt(b.x*b.x+b.y*b.y+b.z*b.z)));
return acos(arg)*180.0/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
stringstream str2ss(string s){
stringstream sstr;
sstr.str(""); sstr.clear(); sstr.str(s);
return sstr;
}
//测试打开输入文件 如果成功则返回0并输出信息 否则返回1
int open_infile(ifstream &infile, const char* filename){
infile.open(filename);
if (!infile){
cerr << BOLDRED << "error ==> " << RESET << "file not found: " << filename << endl;
return -1;
}
return 0;
}
//测试打开输出文件 如果成功则返回0并输出信息 否则返回1
int open_outfile(ofstream &outfile, const char* filename){
outfile.open(filename);
if (!outfile){
cerr << BOLDRED << "error ==> " << RESET << "fail to create the file: " << filename << endl;
return -1;
}
return 0;
}
/*
y
|
|
v3-------v4
| |
| |
| |
v1-------v2-------> x
*/
double cube_interpolate(double x1, double x2, double y1, double y2,
double x, double y, double v1, double v2, double v3, double v4)
{
double A = (x - x1) / (x2 - x1);
double B = (y - y1) / (y2 - y1);
return A * B * (v4 - v3 - v2 + v1) + B * (v3 - v1) + A * (v2 - v1) + v1;
}

181
src/head_func.h Normal file
View File

@@ -0,0 +1,181 @@
#ifndef _HEAD_FUNC_H
#define _HEAD_FUNC_H
#include "iostream"
#include "fstream"
#include "sstream"
#include "string.h"
#include "cmath"
#include "iomanip"
#include "stdio.h"
#include "stdlib.h"
#include "unistd.h"
#include "vector"
#include "map"
#include "algorithm"
#include "ctime"
#include "omp.h"
#include "random"
using namespace std;
//数学常量
#define BDL_MAX 1e+30 ///< 定义变量的最大值
#define BDL_MIN -1e+30 ///< 定义变量的最小值
#define ZERO 1e-20 ///< 定义零值
//物理常量
#define Pi (4.0*atan(1.0)) ///< 圆周率
#define G0 6.67408e-3 ///< 万有引力常数。注意这里的量级本来应该是e-11，考虑到单位转换，取维度单位为m，密度单位为g/cm^3，乘以G0则重力单位即为mGal
#define T0 5e+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; ///< 整型到整形的映射
/**
* @brief 点结构体
*
* 直角坐标系下的一个实点，包含该点的坐标位置与点的编号。
*/
struct cpoint{
int id = -1; ///< 点的索引值。初始值为-1，一般从0开始标号。
double x = BDL_MAX; ///< 点的x坐标值。
double y = BDL_MAX; ///< 点的y坐标值。
double z = BDL_MAX; ///< 点的z坐标值。
};
/**
*@brief cpoint结构体的向量。
*/
typedef vector<cpoint> cpointArray;
/**
* @brief 观测点结构体
*
* 直角坐标系下的一个观测点。该结构体由cpoint继承而来，增加了观测点的数值val与不确定度dev。
*/
struct obspoint : cpoint{
double val = BDL_MAX; ///< 标量观测值，初始值为BDL_MAX。
double dev = BDL_MAX; ///< 观测值的不确定度，初始值为BDL_MAX。
};
/**
* @brief obspoint结构体的向量。
*/
typedef vector<obspoint> obspointArray;
/**
* @brief 方块结构体
*
* 直角坐标系下表示一个方块的结构体。
*/
struct block{
int id = -1; ///< 块体的索引。初始值为-1，一般从0开始编号。
int phy_group = -1; ///< 块体所属的物理组编号。初始值为-1。
int vert[4] = {-1,-1,-1,-1}; ///< 块体的四个角点的索引值。初始值均为-1。
double xmin = BDL_MAX; ///< 块体在x坐标轴的最小值。
double xmax = BDL_MIN; ///< 块体在x坐标轴的最大值。
double ymin = BDL_MAX; ///< 块体在y坐标轴的最小值。
double ymax = BDL_MIN; ///< 块体在y坐标轴的最大值。
double zmin = BDL_MAX; ///< 块体在z坐标轴的最小值。
double zmax = BDL_MIN; ///< 块体在z坐标轴的最大值。
};
/**
* @brief block结构体的向量
*/
typedef vector<block> blockArray;
/*************************矢量函数********************************/
/**
* @brief 矢量减法
*
* @param[in] <unnamed> 输入向量
* @param[in] <unnamed> 输入向量
*
* @return 两个向量的差值。
*/
cpoint operator -(cpoint,cpoint);
/**
* @brief 矢量数乘法
*
* @param[in] <unnamed> 输入标量值
* @param[in] <unnamed> 输入向量
*
* @return 标量与向量的乘积
*/
cpoint operator *(double,cpoint);
/**
* @brief 求一个矢量的模
*
* @param[in] <unnamed> 输入向量
*
* @return 向量的模
*/
double modCpoint(cpoint);
/**
* @brief 求两个向量的夹角
*
* @param[in] <unnamed> 输入向量
* @param[in] <unnamed> 输入向量
*
* @return 输入向量的夹角
*/
double angCpoint(cpoint,cpoint);
/*************************全局函数********************************/
/**
* @brief 返回反正切角
*
* @param[in] <unnamed> 一个double类型的浮点值
*
* @return 一个double类型的浮点值
*/
double arctg(double);
/**
* @brief 将string转换为stringstream
*
* @param[in] <unnamed> 输入字符串
*
* @return CString类型字符串
*/
stringstream str2ss(string);
/**
* @brief 打开输入文件并返回执行情况
*
* @param <unnamed> ifstream对象
* @param <unnamed> 文件名
*
* @return 执行成功返回0，不成功则返回-1并输出错误信息。
*/
int open_infile(ifstream&, const char*);
/**
* @brief 打开输出文件并返回执行情况
*
* @param <unnamed> ofstream对象
* @param <unnamed> 文件名
*
* @return 执行成功返回0，不成功则返回-1并输出错误信息。
*/
int open_outfile(ofstream&, const char*);
/**
* @brief 规则网格内插值
*
* @param[in] <unnamed> 输入参数
*
* @return 给定点位的插值
*/
double cube_interpolate(double, double, double, double, double, double, double, double, double, double); //四方块内的双线性插值
#endif

150
src/init_mod_blocks.cpp Normal file
View File

@@ -0,0 +1,150 @@
#include "sgm_esl.h"
int SGM_ESL::InitModBlocks(char* filename,int region_phy_group,double upper_depth,double lower_depth){
int temp_int,temp_id,ele_type,attri_num;
double temp_val;
_1dArray temp_element;
string temp_str;
stringstream temp_ss;
int2intMap vert_id_map;
ifstream mshin;
if (open_infile(mshin,filename)) return -1; //检查并打开模型文件
while(getline(mshin,temp_str)){
//读入模型空间顶点集 msh文件版本为2.2
if (temp_str == "$Nodes"){
getline(mshin,temp_str);
temp_ss = str2ss(temp_str);
temp_ss >> mod_vertice_num_; //第一个数为顶点的个数
mod_vertice_.resize(mod_vertice_num_); //开辟空间
for (int i = 0; i < mod_vertice_num_; i++){
getline(mshin,temp_str);
temp_ss = str2ss(temp_str);
temp_ss >> mod_vertice_[i].id >> mod_vertice_[i].x >> mod_vertice_[i].y >> mod_vertice_[i].z;
//记录顶点索引对应排序位置用于清理块体顶点索引
vert_id_map[mod_vertice_[i].id] = i;
//重置顶点索引
mod_vertice_[i].id = i;
}
}
//读入模型空间单元体
else if (temp_str == "$Elements"){
getline(mshin,temp_str);
temp_ss = str2ss(temp_str);
temp_ss >> mod_blocks_num_;
mod_blocks_.resize(mod_blocks_num_);
for (int i = 0; i < mod_blocks_num_; i++){
getline(mshin,temp_str);
temp_ss = str2ss(temp_str);
temp_ss >> mod_blocks_[i].id >> ele_type;
if (ele_type != 3){
cerr << BOLDRED << "error ==> " << RESET << "wrong model element type, quit program!" << endl;
return 1;
}
else{
//读入模型单元物理组信息用于分组 注意这里attri_num应该永远为1
temp_ss >> attri_num;
if (attri_num != 1){
cerr << BOLDRED << "error ==> " << RESET << "wrong model element attributes, quit program!" << endl;
return 1;
}
else{
temp_ss >> mod_blocks_[i].phy_group;
for (int j = 0; j < 4; j++){
temp_ss >> mod_blocks_[i].vert[j];
}
}
}
}
}
else continue; //不能识别的单元都被忽略了
}
mshin.close();
//清理模型块体顶点索引排序
for (int i = 0; i < mod_blocks_.size(); i++){
for (int j = 0; j < 4; j++){
mod_blocks_[i].vert[j] = vert_id_map[mod_blocks_[i].vert[j]];
}
}
//第二次读入模型文件 初始化模型单元属性
if (open_infile(mshin,filename)) return -1; //检查并打开模型文件
while(getline(mshin,temp_str)){
if (temp_str == "$ElementData"){
//先读入元素块的名称 添加到数组
for (int i = 0; i < 2; i++)
getline(mshin,temp_str);
mod_elements_names_.push_back(temp_str);
//跳过元素属性前面的值 最后一次为当前元素块的个数
for (int i = 0; i < 6; i++)
getline(mshin,temp_str);
temp_ss = str2ss(temp_str);
temp_ss >> temp_int;
//初始化模型属性
if (!temp_element.empty()) temp_element.clear();
temp_element.resize(mod_blocks_num_,0.0);
for (int i = 0; i < temp_int; i++){
getline(mshin,temp_str);
temp_ss = str2ss(temp_str);
temp_ss >> temp_id >> temp_val; //读入单元体索引与属性值
temp_element[temp_id] = temp_val; //添加单元体属性到数组对应位置
}
mod_elements_.push_back(temp_element);
}
else continue; //不能识别的单元都被忽略了
}
mshin.close();
//模型读入完成 初始化模型单元设置
init_model_.resize(mod_blocks_num_,0.0);
invert_model_.resize(mod_blocks_num_,0.0);
for (int i = 0; i < mod_elements_names_.size(); i++){
if (!strcmp(mod_elements_names_[i].c_str(),init_model_name_)){
init_model_ = mod_elements_[i];
clog << "initial model initialized from file." << endl;
}
else cerr << BOLDYELLOW << "warning ==> " << RESET << "element data not used: " << mod_elements_names_[i] << endl;
}
//计算模型单元尺寸
for (int i = 0; i < mod_blocks_num_; i++){
for (int j = 0; j < 4; j++){
mod_blocks_[i].xmin = MIN(mod_blocks_[i].xmin, mod_vertice_[mod_blocks_[i].vert[j]].x);
mod_blocks_[i].xmax = MAX(mod_blocks_[i].xmax, mod_vertice_[mod_blocks_[i].vert[j]].x);
mod_blocks_[i].ymin = MIN(mod_blocks_[i].ymin, mod_vertice_[mod_blocks_[i].vert[j]].y);
mod_blocks_[i].ymax = MAX(mod_blocks_[i].ymax, mod_vertice_[mod_blocks_[i].vert[j]].y);
}
mod_blocks_[i].zmin = upper_depth;
mod_blocks_[i].zmax = lower_depth;
}
//初始化区域与剩余模型单元列表
regional_model_list_.resize(mod_blocks_num_,0.0);
object_model_list_.resize(mod_blocks_num_,1.0);
if (region_phy_group != -1)
{
for (int i = 0; i < mod_blocks_num_; i++){
if (mod_blocks_[i].phy_group == region_phy_group){
regional_model_list_[i] = 1.0;
object_model_list_[i] = 0.0;
}
}
}
//清理数组
for (int i = 0; i < mod_elements_.size(); i++)
{
mod_elements_[i].clear();
vector <double>().swap(mod_elements_[i]);
}
mod_elements_.clear();
vector < vector<double> >().swap(mod_elements_);
mod_elements_names_.clear();
vector <string>().swap(mod_elements_names_);
temp_element.clear();
vector <double>().swap(temp_element);
return 0;
}

136
src/init_obs_points.cpp Normal file
View File

@@ -0,0 +1,136 @@
#include "sgm_esl.h"
/**
* @brief 初始化观测点数组
*
* 注意我们的模型空间为z轴向下的右手直角坐标系，因此y轴的正方向为东向，而x轴的正方向为北向。这与通常二维数据的x与y轴朝向相反。
*
* @param para_char 观测数据文件与参数
*
* @return 返回执行情况。正确返回0，否则返回-1
*/
int SGM_ESL::InitObsPoints(char* para_char){
bool if_interpolate = false;
char filename[1024];
double obs_xmin, obs_xmax, obs_ymin, obs_ymax, obs_dx, obs_dy;
int obs_xnum, obs_ynum;
obspoint temp_obs;
string temp_str;
//解析文件名+r<xmin>/<xmax>/<ymin>/<ymax>+i<dx>/dy
if (7 != sscanf(para_char, "%[^+]+r%lf/%lf/%lf/%lf+i%lf/%lf",
filename, &obs_ymin, &obs_ymax, &obs_xmin, &obs_xmax, &obs_dy, &obs_dx)){
strcpy(filename,para_char);
}
else{
clog << BOLDYELLOW << "warning ==> " << RESET << "observations are recognized as a regular grid." << endl;
if_interpolate = true;
obs_ynum = (obs_ymax - obs_ymin)/obs_dy + 1;
obs_xnum = (obs_xmax - obs_xmin)/obs_dx + 1;
}
ifstream infile;
if (open_infile(infile,filename)) return -1;
obspointArray temp_obs_points;
while(getline(infile,temp_str)){
if (*(temp_str.begin()) == '#') continue;
else{
//按每行5个数据解析 初始化为含观测值与不确定度的观测点
if (5 == sscanf(temp_str.c_str(),"%lf %lf %lf %lf %lf",
&temp_obs.y,&temp_obs.x,&temp_obs.z,&temp_obs.val,&temp_obs.dev)){
temp_obs.z *= -1.0;
temp_obs.id = temp_obs_points.size();
temp_obs_points.push_back(temp_obs);
}
else{
cerr << BOLDYELLOW << "ignored ==> " << RESET << "wrong input: " << temp_str << endl;
continue;
}
}
}
infile.close();
// 如果读入的观测点数为空则返回错误
if (temp_obs_points.empty())
{
cerr << BOLDRED << "error ==> " << RESET << "fail to initial observations with the parameter: " << filename << endl;
return -1;
}
// 根据所使用的等效源设置将观测数据插值至每个等效源块体的平面中心位置 以达到适当加密或抽稀观测数据的目的
// 注意这里我们默认输入的观测点数据为规则网数据 数据的排序应为由左下到右上
else if (if_interpolate){
int temp_m, temp_n;
obs_points_.reserve(mod_blocks_num_);
for (int i = 0; i < mod_blocks_num_; i++){
temp_obs.id = obs_points_.size();
temp_obs.y = 0.5 * (mod_blocks_[i].ymin + mod_blocks_[i].ymax);
temp_obs.x = 0.5 * (mod_blocks_[i].xmin + mod_blocks_[i].xmax);
temp_m = floor((temp_obs.y - obs_ymin)/obs_dy);
temp_n = floor((temp_obs.x - obs_xmin)/obs_dx);
if (temp_m >= 0 && temp_n >= 0 && temp_m <= obs_ynum - 2 && temp_n <= obs_xnum - 2){
temp_obs.z = cube_interpolate(obs_ymin + temp_m * obs_dy, obs_ymin + (temp_m + 1) * obs_dy,
obs_xmin + temp_n * obs_dx, obs_xmin + (temp_n + 1) * obs_dx,
temp_obs.y, temp_obs.x,
temp_obs_points[temp_n * obs_ynum + temp_m].z,
temp_obs_points[temp_n * obs_ynum + temp_m + 1].z,
temp_obs_points[(temp_n + 1) * obs_ynum + temp_m].z,
temp_obs_points[(temp_n + 1) * obs_ynum + temp_m + 1].z);
temp_obs.val = cube_interpolate(obs_ymin + temp_m * obs_dy, obs_ymin + (temp_m + 1) * obs_dy,
obs_xmin + temp_n * obs_dx, obs_xmin + (temp_n + 1) * obs_dx,
temp_obs.y, temp_obs.x,
temp_obs_points[temp_n * obs_ynum + temp_m].val,
temp_obs_points[temp_n * obs_ynum + temp_m + 1].val,
temp_obs_points[(temp_n + 1) * obs_ynum + temp_m].val,
temp_obs_points[(temp_n + 1) * obs_ynum + temp_m + 1].val);
temp_obs.dev = cube_interpolate(obs_ymin + temp_m * obs_dy, obs_ymin + (temp_m + 1) * obs_dy,
obs_xmin + temp_n * obs_dx, obs_xmin + (temp_n + 1) * obs_dx,
temp_obs.y, temp_obs.x,
temp_obs_points[temp_n * obs_ynum + temp_m].dev,
temp_obs_points[temp_n * obs_ynum + temp_m + 1].dev,
temp_obs_points[(temp_n + 1) * obs_ynum + temp_m].dev,
temp_obs_points[(temp_n + 1) * obs_ynum + temp_m + 1].dev);
obs_points_.push_back(temp_obs);
}
}
obs_points_num_ = obs_points_.size();
//赋值给total_field_ 同时初始化其他数组大小
total_field_.resize(obs_points_num_, 0.0);
regional_field_.resize(obs_points_num_, 0.0);
object_field_.resize(obs_points_num_, 0.0);
regional_regional_field_.resize(obs_points_num_, 0.0);
regional_object_field_.resize(obs_points_num_, 0.0);
object_regional_field_.resize(obs_points_num_, 0.0);
object_object_field_.resize(obs_points_num_, 0.0);
for (int i = 0; i < obs_points_num_; i++){
total_field_[i] = obs_points_[i].val;
}
}
else{
obs_points_ = temp_obs_points;
obs_points_num_ = obs_points_.size();
//赋值给total_field_ 同时初始化其他数组大小
total_field_.resize(obs_points_num_,0.0);
regional_field_.resize(obs_points_num_,0.0);
object_field_.resize(obs_points_num_,0.0);
regional_regional_field_.resize(obs_points_num_,0.0);
regional_object_field_.resize(obs_points_num_,0.0);
object_regional_field_.resize(obs_points_num_,0.0);
object_object_field_.resize(obs_points_num_,0.0);
for (int i = 0; i < obs_points_num_; i++){
total_field_[i] = obs_points_[i].val;
}
}
// 清理临时数组
temp_obs_points.clear();
vector<obspoint>().swap(temp_obs_points);
return 0;
}

11
src/main.cpp Normal file
View File

@@ -0,0 +1,11 @@
#include "sgm_esl.h"
int main(int argc, char* argv[]){
SGM_ESL instance;
if (argc == 0)
{
clog << "usage: GravSepa.ex <para-file>" << endl;
}
else instance.Routine(argv[1]);
return 0;
}

115
src/optimize_cg_precondition.cpp Normal file
View File

@@ -0,0 +1,115 @@
#include "sgm_esl.h"
void SGM_ESL::OptimizeCG_Precondition(_1dArray obs_array){
int i,j;
double data_object_func, dTAd, zTr, z1Tr1, ak, beta_k;
_1dArray r_k, r_k1, z_k, z_k1, d_k;
_1dArray Adk, Adk_part;
_1dArray predict_field;
r_k.resize(mod_blocks_num_); r_k1.resize(mod_blocks_num_);
z_k.resize(mod_blocks_num_); z_k1.resize(mod_blocks_num_);
d_k.resize(mod_blocks_num_);
Adk.resize(mod_blocks_num_); Adk_part.resize(obs_points_num_);
predict_field.resize(obs_points_num_);
//初始化反演密度向量
for (i = 0; i < mod_blocks_num_; i++)
invert_model_[i] = init_model_[i];
#pragma omp parallel for private(i,j) schedule(guided)
for (i = 0; i < obs_points_num_ ; i++){
Adk_part[i] = 0.0;
for (j = 0; j < mod_blocks_num_; j++){
Adk_part[i] += GM_kernel_[i][j]*invert_model_[j];
}
}
#pragma omp parallel for private(i,j) schedule(guided)
for (j = 0; j < mod_blocks_num_; j++){
Adk[j] = 0.0;
for (i = 0; i < obs_points_num_ ; i++){
Adk[j] += GM_kernel_[i][j]*WdTWd[i]*Adk_part[i];
}
}
#pragma omp parallel for private(i) schedule(guided)
for (j = 0; j < mod_blocks_num_; j++)
r_k[j] = PartB[j] - Adk[j];
#pragma omp parallel for private(i) schedule(guided)
for (j = 0; j < mod_blocks_num_; j++){
z_k[j] = Wp[j]*r_k[j];
d_k[j] = z_k[j];
}
for (int t = 0; t < iter_times_; t++){
//计算数据拟合差函数
#pragma omp parallel for private(i,j) schedule(guided)
for (i = 0; i < obs_points_num_; i++){
predict_field[i] = 0.0;
for (j = 0; j < mod_blocks_num_; j++){
predict_field[i] += GM_kernel_[i][j] * invert_model_[j];
}
}
data_object_func = 0.0;
for (i = 0; i < obs_points_num_; i++){
data_object_func += pow((obs_array[i] - predict_field[i])/obs_points_[i].dev,2);
}
data_object_func /= obs_points_num_;
if (data_object_func <= 1.0){
cout << "Iteration limit reached!\tconvergence value = " << data_object_func << endl;
cout << "===================================" << endl;
break;
}
else{
cout << "Iteration times: " << t << "\tconvergence value = " << data_object_func << endl;
MOVEUP(1);
cout << CLEARLINE;
}
//计算r_0
#pragma omp parallel for private(i,j) schedule(guided)
for (i = 0; i < obs_points_num_ ; i++){
Adk_part[i] = 0.0;
for (j = 0; j < mod_blocks_num_; j++){
Adk_part[i] += GM_kernel_[i][j]*d_k[j];
}
}
#pragma omp parallel for private(i,j) schedule(guided)
for (j = 0; j < mod_blocks_num_; j++){
Adk[j] = 0.0;
for (i = 0; i < obs_points_num_ ; i++){
Adk[j] += GM_kernel_[i][j]*WdTWd[i]*Adk_part[i];
}
}
zTr = 0.0; dTAd = 0.0;
for (j = 0; j < mod_blocks_num_; j++){
zTr += z_k[j]*r_k[j];
dTAd += d_k[j]*Adk[j];
}
ak = zTr/dTAd;
for (j = 0; j < mod_blocks_num_; j++){
invert_model_[j] += ak*d_k[j];
r_k1[j] = r_k[j] - ak*Adk[j];
z_k1[j] = Wp[j]*r_k1[j];
}
z1Tr1 = 0.0;
for (j = 0; j < mod_blocks_num_; j++)
z1Tr1 += z_k1[j]*r_k1[j];
beta_k = z1Tr1/zTr;
for (j = 0; j < mod_blocks_num_; j++){
d_k[j] = z_k1[j] + beta_k*d_k[j];
r_k[j] = r_k1[j];
z_k[j] = z_k1[j];
}
}
return;
}

25
src/out_mod_blocks.cpp Normal file
View File

@@ -0,0 +1,25 @@
#include "sgm_esl.h"
int SGM_ESL::OutModBlocks(char* filename){
if (!strcmp(filename,"NULL")) return 0;
ofstream outfile;
if (open_outfile(outfile,filename)) return -1;
outfile<<"$MeshFormat"<<endl<<"2.2 0 8"<<endl<<"$EndMeshFormat"<<endl<<"$Nodes"<<endl<<mod_vertice_num_<<endl;
for (int i = 0; i < mod_vertice_num_; i++){
outfile << mod_vertice_[i].id << " " << setprecision(16) << mod_vertice_[i].x << " " << mod_vertice_[i].y << " " << mod_vertice_[i].z << endl;
}
outfile<<"$EndNodes"<<endl;
outfile<<"$Elements"<<endl<<mod_blocks_num_<<endl;
for (int i = 0; i < mod_blocks_num_; i++){
outfile << mod_blocks_[i].id << " 3 1 " << mod_blocks_[i].phy_group;
for (int j = 0; j < 4; j++){
outfile << " " << mod_blocks_[i].vert[j];
}
outfile << endl;
}
outfile << "$EndElements"<< endl;
outfile.close();
return 0;
}

22
src/out_mod_ele_data.cpp Normal file
View File

@@ -0,0 +1,22 @@
#include "sgm_esl.h"
int SGM_ESL::OutModEleData(char* filename,string data_name,_1dArray data_array){
if (!strcmp(filename,"NULL")) return 0;
ofstream outfile(filename,ios::app);
if (!outfile){
cerr << BOLDRED << "error ==> " << RESET << "fail to create the file: " << filename << endl;
return -1;
}
outfile<<"$ElementData"<<endl;
outfile<<1<<endl<<data_name<<endl<<1<<endl<<0.0<<endl<<3<<endl<<0<<endl<<1<<endl<<data_array.size()<<endl;
for (int i = 0; i < data_array.size(); i++)
{
//注意这里我们没有记录稀疏模型的id 所以不适合输出稀疏模型
outfile << i << " " << setprecision(16) << data_array[i] << endl;
}
outfile<<"$EndElementData"<< endl;
outfile.close();
return 0;
}

31
src/out_obs_points.cpp Normal file
View File

@@ -0,0 +1,31 @@
#include "sgm_esl.h"
int SGM_ESL::OutObsPoints(const char* filename,_1dArray input_data){
ofstream outfile;
if (open_outfile(outfile,filename)) return -1;
outfile << "# y-northing(m) x-easting(m) ele(m) obs-val(mGal|nT|Eo) obs-dev(mGal|nT|Eo) reproduced-data(mGal|nT|Eo)" << endl;
for (int i = 0; i < obs_points_num_; i++)
{
outfile << setprecision(16) << obs_points_[i].y << " " << obs_points_[i].x << " "
<< -1.0*obs_points_[i].z << " " << obs_points_[i].val << " " << obs_points_[i].dev << " " << input_data[i] << endl;
}
outfile.close();
return 0;
}
int SGM_ESL::OutObsPoints(const char* filename,_1dArray reg_data,_1dArray loc_data,_1dArray res_data){
ofstream outfile;
if (open_outfile(outfile,filename)) return -1;
outfile << "# y-northing(m) x-easting(m) ele(m) obs-val(mGal|nT|Eo) \
obs-dev(mGal|nT|Eo) total-recovered(mGal|nT|Eo) regional-recovered(mGal|nT|Eo) \
local-recovered(mGal|nT|Eo) residual-recovered(mGal|nT|Eo) difference (mGal|nT|Eo)" << endl;
for (int i = 0; i < obs_points_num_; i++)
{
outfile << setprecision(16) << obs_points_[i].y << " " << obs_points_[i].x << " "
<< -1.0*obs_points_[i].z << " " << obs_points_[i].val << " " << obs_points_[i].dev << " "
<< reg_data[i]+loc_data[i] << " " << reg_data[i] << " " << loc_data[i] << " " << res_data[i] << " "
<< reg_data[i]+loc_data[i] - obs_points_[i].val << endl;
}
outfile.close();
return 0;
}

114
src/progressBar.cpp Normal file
View File

@@ -0,0 +1,114 @@
//#ifdef _WINDOWS
//#include <windows.h>
//#else
//#include <sys/ioctl.h>
//#endif
#include "progressBar.h"
ProgressBar::ProgressBar() {}
ProgressBar::ProgressBar(unsigned long n_, const char* description_, std::ostream& out_){
n = n_;
frequency_update = n_;
description = description_;
out = &out_;
unit_bar = "\u2588";
unit_space = "-";
desc_width = std::strlen(description); // character width of description field
}
void ProgressBar::SetFrequencyUpdate(unsigned long frequency_update_){
if(frequency_update_ > n){
frequency_update = n; // prevents crash if freq_updates_ > n_
}
else{
frequency_update = frequency_update_;
}
}
void ProgressBar::SetStyle(const char* unit_bar_, const char* unit_space_){
unit_bar = unit_bar_;
unit_space = unit_space_;
}
int ProgressBar::GetConsoleWidth(){
int width;
#ifdef _WINDOWS
CONSOLE_SCREEN_BUFFER_INFO csbi;
GetConsoleScreenBufferInfo(GetStdHandle(STD_OUTPUT_HANDLE), &csbi);
width = csbi.srWindow.Right - csbi.srWindow.Left;
#else
struct winsize win;
//注意！当我们使用pipe here-doc等通道获取程序参数时无法正确的获取窗口大小 此时我们将使用预定值
if (ioctl(0, TIOCGWINSZ, &win) != -1)
width = win.ws_col;
else width = 100;
#endif
return width;
}
int ProgressBar::GetBarLength(){
// get console width and according adjust the length of the progress bar
int bar_length = static_cast<int>((GetConsoleWidth() - desc_width - CHARACTER_WIDTH_PERCENTAGE) / 2.);
return bar_length;
}
void ProgressBar::ClearBarField(){
for(int i=0;i<GetConsoleWidth();++i){
*out << " ";
}
*out << "\r" << std::flush;
}
void ProgressBar::Progressed(unsigned long idx_)
{
try{
if(idx_ > n) throw idx_;
// determines whether to update the progress bar from frequency_update
if ((idx_ != n-1) && ((idx_+1) % (n/frequency_update) != 0)) return;
// calculate the size of the progress bar
int bar_size = GetBarLength();
// calculate percentage of progress
double progress_percent = idx_* TOTAL_PERCENTAGE/(n-1);
// calculate the percentage value of a unit bar
double percent_per_unit_bar = TOTAL_PERCENTAGE/bar_size;
// display progress bar
*out << " " << description << " |";
for(int bar_length=0;bar_length<=bar_size-1;++bar_length){
if(bar_length*percent_per_unit_bar<progress_percent){
*out << unit_bar;
}
else{
*out << unit_space;
}
}
if(idx_ == n-1)
*out << "|" << std::setw(CHARACTER_WIDTH_PERCENTAGE + 1) << std::setprecision(1) << std::fixed << progress_percent << "%\r" << std::flush << std::endl;
else *out << "|" << std::setw(CHARACTER_WIDTH_PERCENTAGE + 1) << std::setprecision(1) << std::fixed << progress_percent << "%\r" << std::flush;
}
catch(unsigned long e){
ClearBarField();
std::cerr << "PROGRESS_BAR_EXCEPTION: _idx (" << e << ") went out of bounds, greater than n (" << n << ")." << std::endl << std::flush;
}
}

66
src/progressBar.h Normal file
View File

@@ -0,0 +1,66 @@
#ifndef _PROGRESS_BAR_
#define _PROGRESS_BAR_
#include <sys/ioctl.h>
#include <iostream>
#include <iomanip>
#include <cstring>
#include <thread>
#include <chrono>
#define TOTAL_PERCENTAGE 100.0 //!< 设置总的百分比为100%。
#define CHARACTER_WIDTH_PERCENTAGE 4 //!< 设置每显示一个字符所占的百分比为5%。
/**
* @brief 终端显示进度条的类。
* @warning 这个类是一个在终端显示进度的通用类型，与本程序无直接关系。
*/
class ProgressBar
{
public:
ProgressBar(); ///< 类默认的构造函数。
/**
* @brief 类默认的构造函数。
*
* @param[in] n_ 整形，进度条循环的最大值。
* @param[in] description_ 字符串，进度条的名称。
* @param out_ 输出至标准错误输出。
*/
ProgressBar(unsigned long n_, const char *description_="", std::ostream& out_=std::cerr);
/**
* @brief 设置进度条刷新频率。
*
* @param[in] frequency_update_ 整形，刷新的频率
*/
void SetFrequencyUpdate(unsigned long frequency_update_);
/**
* @brief 设置进度条的风格。
*
* @param[in] unit_bar_ 字符串，已达到进度的符号。
* @param[in] unit_space_ 字符串，未达到进度的符号。
*/
void SetStyle(const char* unit_bar_, const char* unit_space_);
/**
* @brief 更新当前已达到的进度值。
*
* @param[in] idx_ 整形，当前进度条循环的值。
*/
void Progressed(unsigned long idx_);
private:
unsigned long n;
unsigned int desc_width;
unsigned long frequency_update;
std::ostream* out;
const char *description;
const char *unit_bar;
const char *unit_space;
void ClearBarField();
int GetConsoleWidth();
int GetBarLength();
};
#endif

239
src/routine.cpp Normal file
View File

@@ -0,0 +1,239 @@
#include "sgm_esl.h"
int SGM_ESL::Routine(char* parafile){
double object_func, delta_grav;
string temp_str;
//初始化char数组
char run_type[1024], data_type[1024], msh_name[1024], msh_out_name[1024], obs_name[1024], obsout_name[1024];
strcpy(run_type,"NULL");
strcpy(data_type,"NULL");
strcpy(msh_name,"NULL");
strcpy(msh_out_name,"NULL");
strcpy(obs_name,"NULL");
strcpy(obsout_name,"NULL");
strcpy(init_model_name_,"NULL");
//初始化参数
int regional_tag = -1;
double layer_depth = 50, layer_Depth = 100;
alpha_ = 1.0;
beta_ = 1.0;
iter_times_ = 300;
I0_ = loc_I_ = reg_I_ = 90;
D0_ = loc_D_ = reg_D_ = 0;
ifstream infile;
if (open_infile(infile,parafile)) return -1;
while(getline(infile,temp_str)){
if (*(temp_str.begin()) == '#') continue;
else{
if (1 == sscanf(temp_str.c_str(),"run-type=%[^\n]",run_type))
clog << "run type:\t" << run_type << endl;
else if (1 == sscanf(temp_str.c_str(),"data-type=%[^\n]",data_type))
clog << "data type:\t" << data_type << endl;
else if (1==sscanf(temp_str.c_str(),"mshin-file=%[^\n]",msh_name))
clog << "mshin file:\t" << msh_name << endl;
else if (1==sscanf(temp_str.c_str(),"mshout-file=%[^\n]",msh_out_name))
clog << "mshout file:\t" << msh_out_name << endl;
else if (1==sscanf(temp_str.c_str(),"obsin-file=%[^\n]",obs_name))
clog << "obsin file:\t" << obs_name << endl;
else if (1==sscanf(temp_str.c_str(),"obsout-file=%[^\n]",obsout_name))
clog << "obsout file:\t" << obsout_name << endl;
else if (1==sscanf(temp_str.c_str(),"initial-model=%[^\n]",init_model_name_))
clog << "initial model:\t" << init_model_name_ << endl;
else if (1==sscanf(temp_str.c_str(),"depth-weighting=%lf",&beta_))
clog << "depth weighting:\t" << beta_ << endl;
else if (1==sscanf(temp_str.c_str(),"volume-weighting=%lf",&alpha_))
clog << "volume weighting:\t" << alpha_ << endl;
else if (1==sscanf(temp_str.c_str(),"regional-tag=%d",&regional_tag))
clog << "regional tag:\t" << regional_tag << endl;
else if (2==sscanf(temp_str.c_str(),"layer-depth=%lf/%lf",&layer_depth,&layer_Depth))
clog << "layer depth:\t" << layer_depth << "-->" << layer_Depth << endl;
else if (6==sscanf(temp_str.c_str(),"magnetic-para=%lf/%lf/%lf/%lf/%lf/%lf",&I0_,&D0_,&loc_I_,&loc_D_,&reg_I_,&reg_D_))
clog << "incident/decline angle of geomagnetic field:\t" << I0_ << "/" << D0_ << endl
<< "incident/decline angle of magnetization (local):\t" << loc_I_ << "/" << loc_D_ << endl
<< "incident/decline angle of magnetization (regional):\t" << reg_I_ << "/" << reg_D_ << endl;
else if (1==sscanf(temp_str.c_str(),"iteration-para=%d",&iter_times_))
clog << "iteration times: " << iter_times_ << endl;
else
{
cerr << BOLDRED << "error ==> " << RESET << "undefined parameter: " << temp_str << endl;
return -1;
}
}
}
infile.close();
clog << "================================" << endl;
if (InitModBlocks(msh_name,regional_tag,layer_depth,layer_Depth)) return -1;
if (InitObsPoints(obs_name)) return -1;
if (!strcmp(run_type,"inversion")){
//计算核函数
if (!strcmp(data_type,"gravity")) CalGKernel();
else if (!strcmp(data_type,"deltaT")) CalDeltaTKernel();
else{
cerr << BOLDRED << "error ==> " << RESET << "unknown data type: " << data_type << endl;
return -1;
}
//计算加权矩阵
CalWvWz();
CalWp(1.0);
//首先按1.0计算数据拟合差矩阵
CalWd(1.0);
//计算total_feild_对应的PartB
CalPartB(total_field_);
//反演total_field_
OptimizeCG_Precondition(total_field_);
OutModBlocks(msh_out_name);
OutModEleData(msh_out_name,"\"equivalent source\"",invert_model_);
//输出分离的局部与区域场
CalPartialModel(regional_field_);
if (OutObsPoints(obsout_name, regional_field_)) return -1;
}
else if (!strcmp(run_type,"separate")){
//计算核函数
if (!strcmp(data_type,"gravity")) CalGKernel();
else if (!strcmp(data_type,"deltaT")) CalDeltaTKernel();
else{
cerr << BOLDRED << "error ==> " << RESET << "unknown data type: " << data_type << endl;
return -1;
}
//计算加权矩阵
CalWvWz();
CalWp(1.0);
//首先按1.0计算数据拟合差矩阵
CalWd(1.0);
//计算total_feild_对应的PartB
CalPartB(total_field_);
//反演total_field_
OptimizeCG_Precondition(total_field_);
//正演区域和局部异常
CalPartialModel(regional_field_,regional_model_list_);
CalPartialModel(object_field_,object_model_list_);
//按0.5计算数据拟合差矩阵
CalWd(0.5);
//计算regional_feild_对应的PartB
CalPartB(regional_field_);
//反演regional_feild_
OptimizeCG_Precondition(regional_field_);
//正演regional_feild_区域和局部异常
CalPartialModel(regional_regional_field_,regional_model_list_);
CalPartialModel(regional_object_field_,object_model_list_);
//计算object_feild_对应的PartB
CalPartB(object_field_);
//反演object_feild_
OptimizeCG_Precondition(object_field_);
//正演object_feild_区域和局部异常
CalPartialModel(object_regional_field_,regional_model_list_);
CalPartialModel(object_object_field_,object_model_list_);
// 临时变量 用于收集每次迭代中区域与目标异常的变化
/*
int it_times = 0;
string tmp_str, tmp_str2 = "inter_sep_ret_";
stringstream tmp_ss;
string inter_ret_str;
// 临时代码 用于输出每次迭代中区域与目标异常
tmp_ss.clear(); tmp_ss << it_times; tmp_ss >> tmp_str;
inter_ret_str = tmp_str2 + tmp_str + ".txt";
if (OutObsPoints(inter_ret_str.c_str(),regional_field_,object_field_,regional_object_field_)) return -1;
it_times++;
*/
//计算目标函数
object_func = 0.0;
for (int i = 0; i < obs_points_num_; i++){
object_func += pow((regional_object_field_[i] - object_regional_field_[i])/(sqrt(2)*obs_points_[i].dev),2);
}
object_func /= obs_points_num_;
while(object_func > 1.0){
clog << "Objective function = " << object_func << endl;
//合并新的局部与区域场
for (int i = 0; i < obs_points_num_; i++){
regional_field_[i] = regional_regional_field_[i] + object_regional_field_[i];
object_field_[i] = regional_object_field_[i] + object_object_field_[i];
}
//计算regional_feild_对应的PartB
CalPartB(regional_field_);
//反演regional_feild_
OptimizeCG_Precondition(regional_field_);
//正演regional_feild_区域和局部异常
CalPartialModel(regional_regional_field_,regional_model_list_);
CalPartialModel(regional_object_field_,object_model_list_);
//计算object_feild_对应的PartB
CalPartB(object_field_);
//反演object_feild_
OptimizeCG_Precondition(object_field_);
//正演object_feild_区域和局部异常
CalPartialModel(object_regional_field_,regional_model_list_);
CalPartialModel(object_object_field_,object_model_list_);
// 临时代码 用于输出每次迭代中区域与目标异常
/*
tmp_ss.clear(); tmp_ss << it_times; tmp_ss >> tmp_str;
inter_ret_str = tmp_str2 + tmp_str + ".txt";
if (OutObsPoints(inter_ret_str.c_str(),regional_field_,object_field_,regional_object_field_)) return -1;
it_times++;
*/
//计算目标函数
object_func = 0.0;
for (int i = 0; i < obs_points_num_; i++){
object_func += pow((regional_object_field_[i] - object_regional_field_[i])/(sqrt(2)*obs_points_[i].dev),2);
}
object_func /= obs_points_num_;
}
clog << "Objective function = " << object_func << endl;
double total_engery = 0.0;
double mean_engery = 0.0;
for (int i = 0; i < obs_points_num_; i++){
total_engery += pow(regional_object_field_[i],2);
mean_engery += regional_object_field_[i];
}
clog << "Total energy of the remaining field = " << sqrt(total_engery) << endl;
clog << "RMS of misfit (mGal|nT) = " << sqrt(total_engery/obs_points_num_) << endl;
total_engery = 0.0;
mean_engery /= obs_points_num_;
for (int i = 0; i < obs_points_num_; i++){
total_engery += pow(regional_object_field_[i] - mean_engery,2);
}
clog << "SD of misfit (mGal|nT) = " << sqrt(total_engery/obs_points_num_) << endl;
//输出分离的局部与区域场
if (OutObsPoints(obsout_name,regional_field_,object_field_,regional_object_field_)) return -1;
}
else if (!strcmp(run_type,"RTP")){
//计算核函数
if (!strcmp(data_type,"deltaT")) CalDeltaTKernel();
else{
cerr << BOLDRED << "error ==> " << RESET << "unsupported data type: " << data_type << endl;
return -1;
}
//计算加权矩阵
CalWvWz();
CalWp(1.0);
//首先按1.0计算数据拟合差矩阵
CalWd(1.0);
//计算total_feild_对应的PartB
CalPartB(total_field_);
//反演total_field_
OptimizeCG_Precondition(total_field_);
//重新计算正演核矩阵 垂直极化的情况
I0_ = loc_I_ = reg_I_ = 90.0;
D0_ = loc_D_ = reg_D_ = 0.0;
CalDeltaTKernel();
//输出分离的局部与区域场
CalPartialModel(regional_field_);
if (OutObsPoints(obsout_name, regional_field_)) return -1;
}
else cerr << BOLDRED << "error ==> " << RESET << "unknown running type: " << run_type << endl;
return 0;
}

165
src/sgm_esl.h Normal file
View File

@@ -0,0 +1,165 @@
#ifndef _SGM_ESL_H
#define _SGM_ESL_H
#include "head_func.h"
/**
* @brief 利用等效源对重磁场进行区域场与剩余场的分离
*
* 通过对等效源进行分区，迭代反演并正演不同分区的等效源对重磁场的区域与剩余场的分离。
*/
class SGM_ESL
{
public:
SGM_ESL(){}
~SGM_ESL(){}
/**
* @brief 运行流程函数
*
* @param <unnamed> 参数配置文件
*
* @return 返回运行状态，正确0、错误-1。
*/
int Routine(char*);
/**
* @brief 初始化观测点数据
*
* @warning 规则网格数据如果给定数据范围与间距则函数会将观测点插值至等效源块体的中心位置。
*
* @param <unnamed> 观测点读入参数
*
* 参数形式：<filename>[+r<xmin>/<xmax>/<ymin>/<ymax>+i<dx>/<dy>]
*
* @return 返回运行状态，正确0、错误-1。
*/
int InitObsPoints(char*);
/**
* @brief 输出观测点数据
*
* @param <unnamed> 输出文件名
* @param <unnamed> 需要输出的数据
*
* @return 返回运行状态，正确0、错误-1。
*/
int OutObsPoints(const char*,_1dArray);
/**
* @brief 输出观测点数据
*
* @param <unnamed> 输出文件名
* @param <unnamed> 需要输出的数据
* @param <unnamed> 需要输出的数据
* @param <unnamed> 需要输出的数据
*
* @return 返回运行状态，正确0、错误-1。
*/
int OutObsPoints(const char*,_1dArray,_1dArray,_1dArray);
/**
* @brief 初始化等效源块体
*
* @param <unnamed> 输入文件名
* @param <unnamed> 区域场对应的等效源块体的物理分组代码
* @param <unnamed> 等效源层的上顶面
* @param <unnamed> 等效源层的下底面
*
* @return 返回运行状态，正确0、错误-1。
*/
int InitModBlocks(char*,int,double,double);
/**
* @brief 输出等效源模型块体的框架信息
*
* @param <unnamed> 输出文件名
*
* @return 返回运行状态，正确0、错误-1。
*/
int OutModBlocks(char*);
/**
* @brief 输出等效源模型块体的模型参数信息
*
* @param <unnamed> 输出文件名
* @param[in] <unnamed> 输出的模型参数的名称
* @param. <unnamed>..输出模型参数的数组
*
* @return 返回运行状态，正确0、错误-1。
*/
int OutModEleData(char*,string,_1dArray);
/**
* @brief 计算重力计算核函数
*/
void CalGKernel();
/**
* @brief 计算磁场DeltaT核函数
*/
void CalDeltaTKernel();
/**
* @brief 计算等效源体积与深度加权矩阵
*/
void CalWvWz();
/**
* @brief 计算反演预优矩阵
*
* @param[in] <unnamed> 预优因子，一般为1.0
*/
void CalWp(double);
/**
* @brief 计算模型拟合差矩阵
*
* @param[in] <unnamed> 数据不确定度参数，一般为1.0
*/
void CalWd(double);
/**
* @brief 计算共轭梯度目标矩阵
*
* @param[in] <unnamed> 观测值数组
*/
void CalPartB(_1dArray);
/**
* @brief 计算部分等效源重力数据
*
* @param <unnamed> 输出观测数据数组
*/
void CalPartialModel(_1dArray&);
/**
* @brief 计算部分等效源重力数据
*
* @param <unnamed> 输出观测数据数组
* @param <unnamed> 等效源分组信息数据（可选）
*/
void CalPartialModel(_1dArray&,_1dArray);
/**
* @brief 预优共轭梯度求解函数
*
* @param[in] <unnamed> 观测值数组
*/
void OptimizeCG_Precondition(_1dArray);
private:
//观测点数组
int obs_points_num_;
obspointArray obs_points_;
_1dArray total_field_;
_1dArray regional_field_;
_1dArray object_field_;
_1dArray regional_regional_field_;
_1dArray regional_object_field_;
_1dArray object_regional_field_;
_1dArray object_object_field_;
//模型块体数组
int mod_vertice_num_;
int mod_blocks_num_;
cpointArray mod_vertice_;
blockArray mod_blocks_;
_2dArray mod_elements_;
_1sArray mod_elements_names_;
//模型属性数组
_1dArray init_model_; char init_model_name_[1024];
_1dArray invert_model_;
_1dArray regional_model_list_;
_1dArray object_model_list_;
_2dArray GM_kernel_;
//反演参数
_1dArray WvTWv,WzTWz,WdTWd,Wp,PartB;
int iter_times_;
double alpha_, beta_;
//double D0_, I0_, D_, I_; ///< 地磁参数与磁化参数
// 假设局部与区域磁化参数不一样
double D0_, I0_, loc_D_, loc_I_, reg_D_, reg_I_;
};
#endif //_SGM_ESL_H