update src
This commit is contained in:
yizhangss
parent
3be92e13b4
commit
a7c1477b2e
@ -1,19 +1,13 @@
|
||||
aux_source_directory(. SRC_DIR)
|
||||
aux_source_directory(. DIR_SRC)
|
||||
|
||||
set(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/bin)
|
||||
|
||||
add_executable(gm3d ${SRC_DIR})
|
||||
include_directories(/usr/local/include)
|
||||
|
||||
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --std=c++11 -O2")
|
||||
add_executable(gm3d ${DIR_SRC})
|
||||
|
||||
# 添加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}")
|
||||
target_link_libraries(gm3d PUBLIC OpenMP::OpenMP_CXX)
|
||||
endif()
|
||||
find_library(GCTL_LIBRARY gctl /usr/local/lib)
|
||||
target_link_libraries(gm3d ${GCTL_LIBRARY})
|
||||
|
||||
set(CMAKE_INSTALL_PREFIX /usr/local)
|
||||
# 安装至${CMAKE_INSTALL_PREFIX}/sbin文件夹
|
||||
install(TARGETS gm3d RUNTIME DESTINATION sbin)
|
||||
@ -1,6 +1,6 @@
|
||||
#include "gm3d.h"
|
||||
//我们读入一个界面数据 插值计算每个块体中心位置的值 然后按情况赋值
|
||||
int GM3D::AddInterfaceBlock(modelist para_list){
|
||||
void GM3D::AddInterfaceBlock(modelist para_list){
|
||||
int m, n, xnum, ynum;
|
||||
double xs,xe,xmin,xmax,dx;
|
||||
double ys,ye,ymin,ymax,dy;
|
||||
@ -9,137 +9,134 @@ int GM3D::AddInterfaceBlock(modelist para_list){
|
||||
string temp_str;
|
||||
stringstream temp_ss;
|
||||
|
||||
cpoint temp_topo;
|
||||
cpointArray face_topo;
|
||||
vertex temp_topo;
|
||||
vertexVector face_topo;
|
||||
|
||||
ifstream infile;
|
||||
char filename[1024];
|
||||
if (1 == sscanf(para_list.mod_para,"%s",filename)){
|
||||
if (open_infile(infile,filename)) return -1;
|
||||
gctl::open_infile(infile, para_list.mod_para);
|
||||
|
||||
while (getline(infile,temp_str)){
|
||||
//#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);
|
||||
ymin = MIN(ys,ye); ymax = MAX(ys,ye);
|
||||
dx = fabs(dx); dy = fabs(dy);
|
||||
xnum = round((xmax - xmin)/dx) + 1;
|
||||
ynum = round((ymax - ymin)/dy) + 1;
|
||||
}
|
||||
else continue;
|
||||
}
|
||||
else{
|
||||
temp_ss = str2ss(temp_str);
|
||||
temp_ss >> temp_topo.x >> temp_topo.y >> temp_topo.z;
|
||||
face_topo.push_back(temp_topo);
|
||||
}
|
||||
}
|
||||
infile.close();
|
||||
|
||||
if (!strcmp(para_list.val_type,"replace/top")){
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
m = floor((model_cube_[i].cen.x - xmin)/dx);
|
||||
n = floor((model_cube_[i].cen.y - ymin)/dy);
|
||||
|
||||
temp_topo.z = grid_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
|
||||
face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
|
||||
|
||||
if (model_cube_[i].cen.z >= temp_topo.z){
|
||||
model_block_val_[i] = para_list.mod_value;
|
||||
model_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (!strcmp(para_list.val_type,"replace/bot")){
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
m = floor((model_cube_[i].cen.x - xmin)/dx);
|
||||
n = floor((model_cube_[i].cen.y - ymin)/dy);
|
||||
|
||||
temp_topo.z = grid_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
|
||||
face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
|
||||
|
||||
if (model_cube_[i].cen.z <= temp_topo.z){
|
||||
model_block_val_[i] = para_list.mod_value;
|
||||
model_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (!strcmp(para_list.val_type,"add/top")){
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
m = floor((model_cube_[i].cen.x - xmin)/dx);
|
||||
n = floor((model_cube_[i].cen.y - ymin)/dy);
|
||||
|
||||
temp_topo.z = grid_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
|
||||
face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
|
||||
|
||||
if (model_cube_[i].cen.z >= temp_topo.z){
|
||||
if (model_block_val_[i] == BDL_MAX)
|
||||
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
|
||||
else
|
||||
model_block_val_[i] += para_list.mod_value;
|
||||
model_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (!strcmp(para_list.val_type,"add/bot")){
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
m = floor((model_cube_[i].cen.x - xmin)/dx);
|
||||
n = floor((model_cube_[i].cen.y - ymin)/dy);
|
||||
|
||||
temp_topo.z = grid_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
|
||||
face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
|
||||
|
||||
if (model_cube_[i].cen.z <= temp_topo.z){
|
||||
if (model_block_val_[i] == BDL_MAX)
|
||||
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
|
||||
else
|
||||
model_block_val_[i] += para_list.mod_value;
|
||||
model_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (!strcmp(para_list.val_type,"erase/top")){
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
m = floor((model_cube_[i].cen.x - xmin)/dx);
|
||||
n = floor((model_cube_[i].cen.y - ymin)/dy);
|
||||
|
||||
temp_topo.z = grid_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
|
||||
face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
|
||||
|
||||
if (model_cube_[i].cen.z >= temp_topo.z){
|
||||
model_block_val_[i] = BDL_MAX;
|
||||
model_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (!strcmp(para_list.val_type,"erase/bot")){
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
m = floor((model_cube_[i].cen.x - xmin)/dx);
|
||||
n = floor((model_cube_[i].cen.y - ymin)/dy);
|
||||
|
||||
temp_topo.z = grid_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
|
||||
face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
|
||||
|
||||
if (model_cube_[i].cen.z <= temp_topo.z){
|
||||
model_block_val_[i] = BDL_MAX;
|
||||
model_added = true;
|
||||
}
|
||||
while (getline(infile,temp_str))
|
||||
{
|
||||
//#range必须出现在数据之前
|
||||
if (temp_str[0] == '#')
|
||||
{
|
||||
if (6 == sscanf(temp_str.c_str(),"# range=%lf/%lf/%lf/%lf/%lf/%lf",
|
||||
&xs,&dx,&xe,&ys,&dy,&ye)){
|
||||
xmin = GCTL_MIN(xs,xe); xmax = GCTL_MAX(xs,xe);
|
||||
ymin = GCTL_MIN(ys,ye); ymax = GCTL_MAX(ys,ye);
|
||||
dx = fabs(dx); dy = fabs(dy);
|
||||
xnum = round((xmax - xmin)/dx) + 1;
|
||||
ynum = round((ymax - ymin)/dy) + 1;
|
||||
}
|
||||
else continue;
|
||||
}
|
||||
else{
|
||||
cerr << BOLDRED << "error ==> " << RESET << "wrong value type: " << para_list.val_type << " of the model type: " << para_list.mod_type << endl;
|
||||
return -1;
|
||||
gctl::str2ss(temp_str, temp_ss);
|
||||
temp_ss >> temp_topo.x >> temp_topo.y >> temp_topo.z;
|
||||
face_topo.push_back(temp_topo);
|
||||
}
|
||||
}
|
||||
else{
|
||||
cerr << BOLDRED << "error ==> " << RESET << "fail to add blocks with the parameter: " << para_list.mod_para << endl;
|
||||
return -1;
|
||||
}
|
||||
infile.close();
|
||||
|
||||
if (!model_added){
|
||||
cerr << BOLDYELLOW << "warning ==> " << RESET << "no block changed with the parameter: " << para_list.mod_para << endl;
|
||||
return -1;
|
||||
if (para_list.val_type == "replace/top")
|
||||
{
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
m = floor((model_cube_[i].cen.x - xmin)/dx);
|
||||
n = floor((model_cube_[i].cen.y - ymin)/dy);
|
||||
|
||||
temp_topo.z = gctl::rect_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
|
||||
face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
|
||||
|
||||
if (model_cube_[i].cen.z >= temp_topo.z){
|
||||
model_block_val_[i] = para_list.mod_value;
|
||||
model_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
return 0;
|
||||
else if (para_list.val_type == "replace/bot")
|
||||
{
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
m = floor((model_cube_[i].cen.x - xmin)/dx);
|
||||
n = floor((model_cube_[i].cen.y - ymin)/dy);
|
||||
|
||||
temp_topo.z = gctl::rect_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
|
||||
face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
|
||||
|
||||
if (model_cube_[i].cen.z <= temp_topo.z){
|
||||
model_block_val_[i] = para_list.mod_value;
|
||||
model_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (para_list.val_type == "add/top")
|
||||
{
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
m = floor((model_cube_[i].cen.x - xmin)/dx);
|
||||
n = floor((model_cube_[i].cen.y - ymin)/dy);
|
||||
|
||||
temp_topo.z = gctl::rect_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
|
||||
face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
|
||||
|
||||
if (model_cube_[i].cen.z >= temp_topo.z){
|
||||
if (model_block_val_[i] == GCTL_BDL_MAX)
|
||||
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
|
||||
else
|
||||
model_block_val_[i] += para_list.mod_value;
|
||||
model_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (para_list.val_type == "add/bot")
|
||||
{
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
m = floor((model_cube_[i].cen.x - xmin)/dx);
|
||||
n = floor((model_cube_[i].cen.y - ymin)/dy);
|
||||
|
||||
temp_topo.z = gctl::rect_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
|
||||
face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
|
||||
|
||||
if (model_cube_[i].cen.z <= temp_topo.z){
|
||||
if (model_block_val_[i] == GCTL_BDL_MAX)
|
||||
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
|
||||
else
|
||||
model_block_val_[i] += para_list.mod_value;
|
||||
model_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (para_list.val_type == "erase/top")
|
||||
{
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
m = floor((model_cube_[i].cen.x - xmin)/dx);
|
||||
n = floor((model_cube_[i].cen.y - ymin)/dy);
|
||||
|
||||
temp_topo.z = gctl::rect_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
|
||||
face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
|
||||
|
||||
if (model_cube_[i].cen.z >= temp_topo.z){
|
||||
model_block_val_[i] = GCTL_BDL_MAX;
|
||||
model_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (para_list.val_type == "erase/bot")
|
||||
{
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
m = floor((model_cube_[i].cen.x - xmin)/dx);
|
||||
n = floor((model_cube_[i].cen.y - ymin)/dy);
|
||||
|
||||
temp_topo.z = gctl::rect_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
|
||||
face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
|
||||
|
||||
if (model_cube_[i].cen.z <= temp_topo.z){
|
||||
model_block_val_[i] = GCTL_BDL_MAX;
|
||||
model_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else throw GM3D_INVALID_MODEL_VAL_TYPE;
|
||||
|
||||
if (!model_added)
|
||||
throw GM3D_NOMODEL_ADDED;
|
||||
return;
|
||||
}
|
||||
@ -1,45 +1,44 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::AddModels(char* filename){
|
||||
void GM3D::AddModels(std::string modentity_file){
|
||||
string temp_str;
|
||||
stringstream temp_ss;
|
||||
modelist temp_list;
|
||||
|
||||
ifstream infile;
|
||||
if (open_infile(infile,filename)) return -1;
|
||||
while(getline(infile,temp_str)){
|
||||
if (*(temp_str.begin()) == '#') continue;
|
||||
gctl::open_infile(infile, modentity_file);
|
||||
|
||||
while(getline(infile, temp_str)){
|
||||
if (temp_str[0] == '#') continue;
|
||||
else{
|
||||
//按每行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);
|
||||
}
|
||||
else{
|
||||
cerr << BOLDYELLOW << "ignored ==> " << RESET << "wrong input: " << temp_str << endl;
|
||||
continue;
|
||||
}
|
||||
gctl::str2ss(temp_str, temp_ss);
|
||||
temp_ss >> temp_list.mod_type >> temp_list.val_type >> temp_list.mod_value
|
||||
>> temp_list.mod_para;
|
||||
|
||||
if (temp_ss.fail())
|
||||
throw GM3D_MODEL_INIT_FAIL;
|
||||
|
||||
model_list_.push_back(temp_list);
|
||||
}
|
||||
}
|
||||
infile.close();
|
||||
|
||||
for (int i = 0; i < model_list_.size(); i++){
|
||||
if (!strcmp(model_list_[i].mod_type,"regular_block")){
|
||||
for (int i = 0; i < model_list_.size(); i++)
|
||||
{
|
||||
if (model_list_[i].mod_type == "regular_block"){
|
||||
AddRegularBlock(model_list_[i]);
|
||||
}
|
||||
else if (!strcmp(model_list_[i].mod_type,"tilted_block")){
|
||||
else if (model_list_[i].mod_type == "tilted_block"){
|
||||
AddTiltedBlock(model_list_[i]);
|
||||
}
|
||||
else if (!strcmp(model_list_[i].mod_type,"sphere")){
|
||||
else if (model_list_[i].mod_type == "sphere"){
|
||||
AddSphereBlock(model_list_[i]);
|
||||
}
|
||||
else if (!strcmp(model_list_[i].mod_type,"interface")){
|
||||
else if (model_list_[i].mod_type == "interface"){
|
||||
AddInterfaceBlock(model_list_[i]);
|
||||
}
|
||||
else{
|
||||
cerr << BOLDYELLOW << "ignored ==> " << RESET << "unknown model type: " << model_list_[i].mod_type << endl;
|
||||
continue;
|
||||
}
|
||||
else throw GM3D_INVALID_MODEL_TYPE;
|
||||
}
|
||||
return 0;
|
||||
return;
|
||||
}
|
||||
@ -1,62 +1,61 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::AddRegularBlock(modelist para_list){
|
||||
void GM3D::AddRegularBlock(modelist para_list){
|
||||
double xs,xe,xmin,xmax;
|
||||
double ys,ye,ymin,ymax;
|
||||
double zs,ze,zmin,zmax;
|
||||
bool model_added = false;
|
||||
|
||||
if (6 == sscanf(para_list.mod_para,"%lf/%lf/%lf/%lf/%lf/%lf",&xs,&xe,&ys,&ye,&zs,&ze)){
|
||||
xmin = MIN(xs,xe); xmax = MAX(xs,xe);
|
||||
ymin = MIN(ys,ye); ymax = MAX(ys,ye);
|
||||
zmin = MIN(zs,ze); zmax = MAX(zs,ze);
|
||||
std::stringstream tmp_ss;
|
||||
gctl::str2ss(para_list.mod_para, tmp_ss, "/");
|
||||
tmp_ss >> xs >> xe >> ys >> ye >> zs >> ze;
|
||||
|
||||
if (!strcmp(para_list.val_type,"replace")){
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
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){
|
||||
if (tmp_ss.fail())
|
||||
throw GM3D_INVALID_MODEL_PARAMETER;
|
||||
|
||||
xmin = GCTL_MIN(xs,xe); xmax = GCTL_MAX(xs,xe);
|
||||
ymin = GCTL_MIN(ys,ye); ymax = GCTL_MAX(ys,ye);
|
||||
zmin = GCTL_MIN(zs,ze); zmax = GCTL_MAX(zs,ze);
|
||||
|
||||
if (para_list.val_type == "replace")
|
||||
{
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
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_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (para_list.val_type == "add")
|
||||
{
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
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){
|
||||
if (model_block_val_[i] == GCTL_BDL_MAX)
|
||||
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
|
||||
model_added = true;
|
||||
}
|
||||
else
|
||||
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++){
|
||||
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){
|
||||
if (model_block_val_[i] == BDL_MAX)
|
||||
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
|
||||
else
|
||||
model_block_val_[i] += para_list.mod_value;
|
||||
model_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (!strcmp(para_list.val_type,"erase")){
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
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_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else{
|
||||
cerr << BOLDRED << "error ==> " << RESET << "wrong value type: " << para_list.val_type << " of the model type: " << para_list.mod_type << endl;
|
||||
return -1;
|
||||
}
|
||||
}
|
||||
else{
|
||||
cerr << BOLDRED << "error ==> " << RESET << "fail to add blocks with the parameter: " << para_list.mod_para << endl;
|
||||
return -1;
|
||||
else if (para_list.val_type == "erase")
|
||||
{
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
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] = GCTL_BDL_MAX; //注意重复赋值的块体会覆盖
|
||||
model_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else throw GM3D_INVALID_MODEL_VAL_TYPE;
|
||||
|
||||
if (!model_added){
|
||||
cerr << BOLDYELLOW << "warning ==> " << RESET << "no block changed with the parameter: " << para_list.mod_para << endl;
|
||||
return -1;
|
||||
}
|
||||
return 0;
|
||||
if (!model_added)
|
||||
throw GM3D_NOMODEL_ADDED;
|
||||
return;
|
||||
}
|
||||
@ -1,82 +1,81 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::AddSphereBlock(modelist para_list){
|
||||
void GM3D::AddSphereBlock(modelist para_list){
|
||||
double xc,yc,zc,rad_x,rad_y,rad_z;
|
||||
double dist,rad_limit;
|
||||
double rel_x,rel_y,rel_z,theta,phi;
|
||||
bool model_added = false;
|
||||
|
||||
if (6 == sscanf(para_list.mod_para,"%lf/%lf/%lf/%lf/%lf/%lf",&xc,&yc,&zc,&rad_x,&rad_y,&rad_z)){
|
||||
if (!strcmp(para_list.val_type,"replace")){
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
rel_x = model_cube_[i].cen.x - xc;
|
||||
rel_y = model_cube_[i].cen.y - yc;
|
||||
rel_z = model_cube_[i].cen.z - zc;
|
||||
dist = sqrt(rel_x*rel_x + rel_y*rel_y + rel_z*rel_z);
|
||||
std::stringstream tmp_ss;
|
||||
gctl::str2ss(para_list.mod_para, tmp_ss, "/");
|
||||
tmp_ss >> xc >> yc >> zc >> rad_x>> rad_y >> rad_z;
|
||||
|
||||
theta = acos(rel_z/dist);
|
||||
phi = atan2(rel_y,rel_x);
|
||||
if (tmp_ss.fail())
|
||||
throw GM3D_INVALID_MODEL_PARAMETER;
|
||||
|
||||
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 (para_list.val_type == "replace")
|
||||
{
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
rel_x = model_cube_[i].cen.x - xc;
|
||||
rel_y = model_cube_[i].cen.y - yc;
|
||||
rel_z = model_cube_[i].cen.z - zc;
|
||||
dist = sqrt(rel_x*rel_x + rel_y*rel_y + rel_z*rel_z);
|
||||
|
||||
if (dist <= rad_limit){
|
||||
theta = acos(rel_z/dist);
|
||||
phi = atan2(rel_y,rel_x);
|
||||
|
||||
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_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (para_list.val_type == "add")
|
||||
{
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
rel_x = model_cube_[i].cen.x - xc;
|
||||
rel_y = model_cube_[i].cen.y - yc;
|
||||
rel_z = model_cube_[i].cen.z - zc;
|
||||
dist = sqrt(rel_x*rel_x + rel_y*rel_y + rel_z*rel_z);
|
||||
|
||||
theta = acos(rel_z/dist);
|
||||
phi = atan2(rel_y,rel_x);
|
||||
|
||||
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){
|
||||
if (model_block_val_[i] == GCTL_BDL_MAX)
|
||||
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
|
||||
model_added = true;
|
||||
}
|
||||
else
|
||||
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++){
|
||||
rel_x = model_cube_[i].cen.x - xc;
|
||||
rel_y = model_cube_[i].cen.y - yc;
|
||||
rel_z = model_cube_[i].cen.z - zc;
|
||||
dist = sqrt(rel_x*rel_x + rel_y*rel_y + rel_z*rel_z);
|
||||
}
|
||||
else if (para_list.val_type == "erase")
|
||||
{
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
rel_x = model_cube_[i].cen.x - xc;
|
||||
rel_y = model_cube_[i].cen.y - yc;
|
||||
rel_z = model_cube_[i].cen.z - zc;
|
||||
dist = sqrt(rel_x*rel_x + rel_y*rel_y + rel_z*rel_z);
|
||||
|
||||
theta = acos(rel_z/dist);
|
||||
phi = atan2(rel_y,rel_x);
|
||||
theta = acos(rel_z/dist);
|
||||
phi = atan2(rel_y,rel_x);
|
||||
|
||||
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));
|
||||
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){
|
||||
if (model_block_val_[i] == BDL_MAX)
|
||||
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
|
||||
else
|
||||
model_block_val_[i] += para_list.mod_value;
|
||||
model_added = true;
|
||||
}
|
||||
if (dist <= rad_limit){
|
||||
model_block_val_[i] = GCTL_BDL_MAX; //注意重复赋值的块体会覆盖
|
||||
model_added = true;
|
||||
}
|
||||
}
|
||||
else if (!strcmp(para_list.val_type,"erase")){
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
rel_x = model_cube_[i].cen.x - xc;
|
||||
rel_y = model_cube_[i].cen.y - yc;
|
||||
rel_z = model_cube_[i].cen.z - zc;
|
||||
dist = sqrt(rel_x*rel_x + rel_y*rel_y + rel_z*rel_z);
|
||||
|
||||
theta = acos(rel_z/dist);
|
||||
phi = atan2(rel_y,rel_x);
|
||||
|
||||
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_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else{
|
||||
cerr << BOLDRED << "error ==> " << RESET << "wrong value type: " << para_list.val_type << " of the model type: " << para_list.mod_type << endl;
|
||||
return -1;
|
||||
}
|
||||
}
|
||||
else{
|
||||
cerr << BOLDRED << "error ==> " << RESET << "fail to add blocks with the parameter: " << para_list.mod_para << endl;
|
||||
return -1;
|
||||
}
|
||||
else throw GM3D_INVALID_MODEL_VAL_TYPE;
|
||||
|
||||
if (!model_added){
|
||||
cerr << BOLDYELLOW << "warning ==> " << RESET << "no block changed with the parameter: " << para_list.mod_para << endl;
|
||||
return -1;
|
||||
}
|
||||
return 0;
|
||||
if (!model_added)
|
||||
throw GM3D_NOMODEL_ADDED;
|
||||
return;
|
||||
}
|
||||
@ -1,6 +1,6 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::AddTiltedBlock(modelist para_list){
|
||||
void GM3D::AddTiltedBlock(modelist para_list){
|
||||
double xs_1,xe_1,xmin_1,xmax_1;
|
||||
double xs_2,xe_2,xmin_2,xmax_2;
|
||||
double ys_1,ye_1,ymin_1,ymax_1;
|
||||
@ -9,78 +9,77 @@ int GM3D::AddTiltedBlock(modelist para_list){
|
||||
double zs,ze,zmin,zmax;
|
||||
bool model_added = false;
|
||||
|
||||
if (10 == sscanf(para_list.mod_para,"%lf/%lf/%lf/%lf/%lf/%lf/%lf/%lf/%lf/%lf",
|
||||
&xs_1,&xe_1,&ys_1,&ye_1,&zs,&xs_2,&xe_2,&ys_2,&ye_2,&ze)){
|
||||
xmin_1 = MIN(xs_1,xe_1); xmax_1 = MAX(xs_1,xe_1);
|
||||
ymin_1 = MIN(ys_1,ye_1); ymax_1 = MAX(ys_1,ye_1);
|
||||
xmin_2 = MIN(xs_2,xe_2); xmax_2 = MAX(xs_2,xe_2);
|
||||
ymin_2 = MIN(ys_2,ye_2); ymax_2 = MAX(ys_2,ye_2);
|
||||
zmin = MIN(zs,ze); zmax = MAX(zs,ze);
|
||||
std::stringstream tmp_ss;
|
||||
gctl::str2ss(para_list.mod_para, tmp_ss, "/");
|
||||
tmp_ss >> xs_1 >> xe_1 >> ys_1 >> ye_1 >> zs
|
||||
>> xs_2 >> xe_2 >> ys_2 >> ye_2 >> ze;
|
||||
|
||||
if (!strcmp(para_list.val_type,"replace")){
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
//计算当前层的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;
|
||||
layer_ymax = (model_cube_[i].cen.z - zmin)*(ymax_2 - ymax_1)/(zmax - zmin) + ymax_1;
|
||||
if (tmp_ss.fail())
|
||||
throw GM3D_INVALID_MODEL_PARAMETER;
|
||||
|
||||
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){
|
||||
xmin_1 = GCTL_MIN(xs_1,xe_1); xmax_1 = GCTL_MAX(xs_1,xe_1);
|
||||
ymin_1 = GCTL_MIN(ys_1,ye_1); ymax_1 = GCTL_MAX(ys_1,ye_1);
|
||||
xmin_2 = GCTL_MIN(xs_2,xe_2); xmax_2 = GCTL_MAX(xs_2,xe_2);
|
||||
ymin_2 = GCTL_MIN(ys_2,ye_2); ymax_2 = GCTL_MAX(ys_2,ye_2);
|
||||
zmin = GCTL_MIN(zs,ze); zmax = GCTL_MAX(zs,ze);
|
||||
|
||||
if (para_list.val_type == "replace")
|
||||
{
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
//计算当前层的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;
|
||||
layer_ymax = (model_cube_[i].cen.z - zmin)*(ymax_2 - ymax_1)/(zmax - zmin) + ymax_1;
|
||||
|
||||
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_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (para_list.val_type == "add")
|
||||
{
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
//计算当前层的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;
|
||||
layer_ymax = (model_cube_[i].cen.z - zmin)*(ymax_2 - ymax_1)/(zmax - zmin) + ymax_1;
|
||||
|
||||
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){
|
||||
if (model_block_val_[i] == GCTL_BDL_MAX)
|
||||
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
|
||||
model_added = true;
|
||||
}
|
||||
else
|
||||
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范围
|
||||
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;
|
||||
layer_ymax = (model_cube_[i].cen.z - zmin)*(ymax_2 - ymax_1)/(zmax - zmin) + ymax_1;
|
||||
}
|
||||
else if (para_list.val_type == "erase")
|
||||
{
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
//计算当前层的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;
|
||||
layer_ymax = (model_cube_[i].cen.z - zmin)*(ymax_2 - ymax_1)/(zmax - zmin) + ymax_1;
|
||||
|
||||
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){
|
||||
if (model_block_val_[i] == BDL_MAX)
|
||||
model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
|
||||
else
|
||||
model_block_val_[i] += para_list.mod_value;
|
||||
model_added = true;
|
||||
}
|
||||
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] = GCTL_BDL_MAX; //注意重复赋值的块体会覆盖
|
||||
model_added = true;
|
||||
}
|
||||
}
|
||||
else if (!strcmp(para_list.val_type,"erase")){
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
//计算当前层的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;
|
||||
layer_ymax = (model_cube_[i].cen.z - zmin)*(ymax_2 - ymax_1)/(zmax - zmin) + ymax_1;
|
||||
}
|
||||
else throw GM3D_INVALID_MODEL_VAL_TYPE;
|
||||
|
||||
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_added = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
else{
|
||||
cerr << BOLDRED << "error ==> " << RESET << "wrong value type: " << para_list.val_type << " of the model type: " << para_list.mod_type << endl;
|
||||
return -1;
|
||||
}
|
||||
}
|
||||
else{
|
||||
cerr << BOLDRED << "error ==> " << RESET << "fail to add blocks with the parameter: " << para_list.mod_para << endl;
|
||||
return -1;
|
||||
}
|
||||
|
||||
if (!model_added){
|
||||
cerr << BOLDYELLOW << "warning ==> " << RESET << "no block changed with the parameter: " << para_list.mod_para << endl;
|
||||
return -1;
|
||||
}
|
||||
return 0;
|
||||
if (!model_added)
|
||||
throw GM3D_NOMODEL_ADDED;
|
||||
return;
|
||||
}
|
||||
@ -1,70 +1,77 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::BuildRegularGrid(char* space_para){
|
||||
cpoint temp_cp;
|
||||
void GM3D::BuildRegularGrid(std::string modspace_para)
|
||||
{
|
||||
vertex temp_cp;
|
||||
cube temp_cu;
|
||||
string temp_id_str;
|
||||
double x,xs,dx,xe,xmin,xmax;
|
||||
double y,ys,dy,ye,ymin,ymax;
|
||||
double z,zs,dz,ze,zmin,zmax;
|
||||
double sign[8][3] = {{-0.5,-0.5,-0.5},{0.5,-0.5,-0.5},{0.5,0.5,-0.5},{-0.5,0.5,-0.5},
|
||||
{-0.5,-0.5,0.5},{0.5,-0.5,0.5},{0.5,0.5,0.5},{-0.5,0.5,0.5}};
|
||||
|
||||
_str2pointMap map_str_point;
|
||||
_str2pointMap::iterator imsp;
|
||||
str2vertMap map_str_point;
|
||||
str2vertMap::iterator imsp;
|
||||
|
||||
if (9 == sscanf(space_para,"%lf/%lf/%lf/%lf/%lf/%lf/%lf/%lf/%lf",
|
||||
&xs,&dx,&xe,&ys,&dy,&ye,&zs,&dz,&ze)){
|
||||
xmin = MIN(xs,xe); xmax = MAX(xs,xe);
|
||||
ymin = MIN(ys,ye); ymax = MAX(ys,ye);
|
||||
zmin = MIN(zs,ze); zmax = MAX(zs,ze);
|
||||
string tmp_str;
|
||||
string tmp_id_str;
|
||||
stringstream tmp_ss;
|
||||
|
||||
temp_cu.dx = fabs(dx); temp_cu.dy = fabs(dy); temp_cu.dz = fabs(dz);
|
||||
y = ys;
|
||||
while(y >= ymin && y <= ymax){
|
||||
x = xs;
|
||||
while(x >= xmin && x <= xmax){
|
||||
z = zs;
|
||||
while(z >= zmin && z <= zmax){
|
||||
//添加denMod
|
||||
temp_cu.cen.id = model_cube_.size();
|
||||
temp_cu.cen.x = x; temp_cu.cen.y = y; temp_cu.cen.z = z;
|
||||
//添加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.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查到当前顶点是否存在,这里需要注意，如果顶点已经存在则只需要将顶点索引置为已存在顶点的索引，不增加顶点计数
|
||||
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;//将新产生的顶点保存到顶点位置映射中
|
||||
}
|
||||
gctl::str2ss(modspace_para, tmp_ss, "/");
|
||||
tmp_ss >> xs >> dx >> xe >> ys >> dy >> ye >> zs >> dz >> ze;
|
||||
|
||||
if (tmp_ss.fail())
|
||||
throw GM3D_INVALID_MESH_PARAMETER;
|
||||
|
||||
xmin = GCTL_MIN(xs,xe); xmax = GCTL_MAX(xs,xe);
|
||||
ymin = GCTL_MIN(ys,ye); ymax = GCTL_MAX(ys,ye);
|
||||
zmin = GCTL_MIN(zs,ze); zmax = GCTL_MAX(zs,ze);
|
||||
temp_cu.dx = fabs(dx);
|
||||
temp_cu.dy = fabs(dy);
|
||||
temp_cu.dz = fabs(dz);
|
||||
|
||||
y = ys;
|
||||
while(y >= ymin && y <= ymax){
|
||||
x = xs;
|
||||
while(x >= xmin && x <= xmax){
|
||||
z = zs;
|
||||
while(z >= zmin && z <= zmax){
|
||||
//添加denMod
|
||||
temp_cu.cen.id = model_cube_.size();
|
||||
temp_cu.cen.x = x; temp_cu.cen.y = y; temp_cu.cen.z = z;
|
||||
//添加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.y = temp_cu.cen.y - sign[i][1]*temp_cu.dy;
|
||||
temp_cp.z = temp_cu.cen.z - sign[i][2]*temp_cu.dz;
|
||||
tmp_id_str = point_id(temp_cp);
|
||||
imsp = map_str_point.find(tmp_id_str);
|
||||
//利用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[tmp_id_str] = temp_cp;//将新产生的顶点保存到顶点位置映射中
|
||||
}
|
||||
model_cube_.push_back(temp_cu);
|
||||
z += dz;
|
||||
}
|
||||
x += dx;
|
||||
model_cube_.push_back(temp_cu);
|
||||
z += dz;
|
||||
}
|
||||
y += dy;
|
||||
x += dx;
|
||||
}
|
||||
y += dy;
|
||||
}
|
||||
|
||||
if (model_cube_.empty()){
|
||||
cerr << BOLDRED << "error ==> " << RESET << "fail to initial model space with the parameter: " << space_para << endl;
|
||||
return -1;
|
||||
}
|
||||
if (model_cube_.empty())
|
||||
throw GM3D_INVALID_MESH_PARAMETER;
|
||||
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_, GCTL_BDL_MAX); //初始化模型块体值为BDL_MAX
|
||||
}
|
||||
return 0;
|
||||
return;
|
||||
}
|
||||
20
src/cpoint_id.cpp
Normal file
20
src/cpoint_id.cpp
Normal file
@ -0,0 +1,20 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
//返回一个cpoint的位置id字符串
|
||||
string GM3D::point_id(vertex c){
|
||||
string vert_id, mid_id;
|
||||
stringstream sstemp;
|
||||
|
||||
sstemp.str(""); sstemp.clear(); sstemp<<setprecision(GCTL_PRECISION)<<c.x;
|
||||
sstemp>>vert_id;
|
||||
|
||||
sstemp.str(""); sstemp.clear(); sstemp<<setprecision(GCTL_PRECISION)<<c.y;
|
||||
sstemp>>mid_id;
|
||||
vert_id = vert_id + " " + mid_id;
|
||||
|
||||
sstemp.str(""); sstemp.clear(); sstemp<<setprecision(GCTL_PRECISION)<<c.z;
|
||||
sstemp>>mid_id;
|
||||
vert_id = vert_id + " " + mid_id;
|
||||
|
||||
return vert_id;
|
||||
}
|
||||
@ -1,329 +0,0 @@
|
||||
#include "disp_help.h"
|
||||
|
||||
void DispHelp::addHeadInfo(string s1,string s2,string s3,string s4)
|
||||
{
|
||||
ex_name = s1; version = s2; descript = s3; author = s4;
|
||||
return;
|
||||
}
|
||||
|
||||
void DispHelp::addUsage(string usg)
|
||||
{
|
||||
usages.push_back(usg);
|
||||
return;
|
||||
}
|
||||
|
||||
void DispHelp::addOption(string msg,string sflag,string lflag)
|
||||
{
|
||||
option tmp_option;
|
||||
tmp_option.message = msg; tmp_option.flag_s = sflag; tmp_option.flag_l = lflag;
|
||||
options.push_back(tmp_option);
|
||||
return;
|
||||
}
|
||||
|
||||
void DispHelp::addOptionSec(string msg,int index)
|
||||
{
|
||||
if (index < 0)
|
||||
{
|
||||
options.back().sec_message.push_back(msg);
|
||||
}
|
||||
else options[index].sec_message.push_back(msg);
|
||||
return;
|
||||
}
|
||||
|
||||
void DispHelp::addExample(string ex)
|
||||
{
|
||||
examples.push_back(ex);
|
||||
return;
|
||||
}
|
||||
|
||||
void DispHelp::changeLayerOut(int left,int right)
|
||||
{
|
||||
front_space = left; back_space = right;
|
||||
return;
|
||||
}
|
||||
|
||||
void DispHelp::show()
|
||||
{
|
||||
int line_length;
|
||||
string segment,full_message;
|
||||
stringstream ss_message;
|
||||
//获取终端窗口的行列数
|
||||
struct winsize w;
|
||||
ioctl(STDOUT_FILENO, TIOCGWINSZ, &w);
|
||||
//显示头信息
|
||||
full_message = ex_name + " " + version + " - " + descript;
|
||||
ss_message.clear(); ss_message.str(full_message);
|
||||
|
||||
line_length = front_space + back_space;
|
||||
while(ss_message >> segment)
|
||||
{
|
||||
if ((line_length+segment.length()+1) <= w.ws_col)
|
||||
{
|
||||
if (line_length == (front_space + back_space))
|
||||
{
|
||||
for (int i = 0; i < front_space; i++) clog << " ";
|
||||
clog << segment << " ";
|
||||
line_length += (segment.length()+1);
|
||||
}
|
||||
else
|
||||
{
|
||||
clog << segment << " ";
|
||||
line_length += (segment.length()+1);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
clog << endl;
|
||||
for (int i = 0; i < front_space+4; i++) clog << " ";
|
||||
clog << segment << " ";
|
||||
line_length = (segment.length()+9+front_space+back_space);
|
||||
}
|
||||
}
|
||||
clog << endl;
|
||||
|
||||
ss_message.clear(); ss_message.str(author);
|
||||
line_length = front_space + back_space;;
|
||||
while(ss_message >> segment)
|
||||
{
|
||||
if ((line_length+segment.length()+1) <= w.ws_col)
|
||||
{
|
||||
if (line_length == (front_space + back_space))
|
||||
{
|
||||
for (int i = 0; i < front_space; i++) clog << " ";
|
||||
clog << "Author: " << segment << " ";
|
||||
line_length += (segment.length()+9);
|
||||
}
|
||||
else
|
||||
{
|
||||
clog << segment << " ";
|
||||
line_length += (segment.length()+1);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
clog << endl;
|
||||
for (int i = 0; i < front_space+4; i++) clog << " ";
|
||||
clog << segment << " ";
|
||||
line_length = (segment.length()+9+front_space+back_space);
|
||||
}
|
||||
}
|
||||
clog << endl;
|
||||
|
||||
if (!usages.empty())
|
||||
{
|
||||
for (int i = 0; i < front_space; i++) clog << " ";
|
||||
clog << "Usage:" << endl;
|
||||
for (int i = 0; i < usages.size(); i++)
|
||||
{
|
||||
ss_message.clear(); ss_message.str(usages[i]);
|
||||
|
||||
line_length = front_space + back_space + 4;
|
||||
while(ss_message >> segment)
|
||||
{
|
||||
if ((line_length+segment.length()+1) <= w.ws_col)
|
||||
{
|
||||
if (line_length == (front_space+back_space+4))
|
||||
{
|
||||
for (int i = 0; i < front_space+4; i++) clog << " ";
|
||||
clog << segment << " ";
|
||||
line_length += (segment.length()+1);
|
||||
}
|
||||
else
|
||||
{
|
||||
clog << segment << " ";
|
||||
line_length += (segment.length()+1);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
clog << endl;
|
||||
for (int i = 0; i < front_space+9; i++) clog << " ";
|
||||
clog << segment << " ";
|
||||
line_length = (segment.length()+10+front_space+back_space);
|
||||
}
|
||||
}
|
||||
clog << endl;
|
||||
}
|
||||
}
|
||||
|
||||
if (!options.empty())
|
||||
{
|
||||
for (int i = 0; i < front_space; i++) clog << " ";
|
||||
clog << "Options:" << endl;
|
||||
for (int i = 0; i < options.size(); i++)
|
||||
{
|
||||
if (options[i].flag_l == "")
|
||||
{
|
||||
full_message = options[i].flag_s+" : "+options[i].message;
|
||||
ss_message.clear(); ss_message.str(full_message);
|
||||
|
||||
line_length = front_space + back_space + 4;
|
||||
while(ss_message >> segment)
|
||||
{
|
||||
if ((line_length+segment.length()+1) <= w.ws_col)
|
||||
{
|
||||
if (line_length == (front_space+back_space+4))
|
||||
{
|
||||
for (int i = 0; i < front_space+4; i++) clog << " ";
|
||||
clog << segment << " ";
|
||||
line_length += (segment.length()+1);
|
||||
}
|
||||
else
|
||||
{
|
||||
clog << segment << " ";
|
||||
line_length += (segment.length()+1);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
clog << endl;
|
||||
for (int i = 0; i < front_space+9; i++) clog << " ";
|
||||
clog << segment << " ";
|
||||
line_length = (segment.length()+10+front_space+back_space);
|
||||
}
|
||||
}
|
||||
clog << endl;
|
||||
|
||||
if (!options[i].sec_message.empty())
|
||||
{
|
||||
for (int j = 0; j < options[i].sec_message.size(); j++)
|
||||
{
|
||||
ss_message.clear(); ss_message.str(options[i].sec_message[j]);
|
||||
|
||||
line_length = front_space + back_space + 9;
|
||||
while(ss_message >> segment)
|
||||
{
|
||||
if ((line_length+segment.length()+1) <= w.ws_col)
|
||||
{
|
||||
if (line_length == (front_space+back_space+9))
|
||||
{
|
||||
for (int i = 0; i < front_space+9; i++) clog << " ";
|
||||
clog << segment << " ";
|
||||
line_length += (segment.length()+1);
|
||||
}
|
||||
else
|
||||
{
|
||||
clog << segment << " ";
|
||||
line_length += (segment.length()+1);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
clog << endl;
|
||||
for (int i = 0; i < front_space+13; i++) clog << " ";
|
||||
clog << segment << " ";
|
||||
line_length = (segment.length()+14+front_space+back_space);
|
||||
}
|
||||
}
|
||||
clog << endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
full_message = options[i].flag_s+" | "+options[i].flag_l+" : "+options[i].message;
|
||||
ss_message.clear(); ss_message.str(full_message);
|
||||
|
||||
line_length = front_space + back_space + 4;
|
||||
while(ss_message >> segment)
|
||||
{
|
||||
if ((line_length+segment.length()+1) <= w.ws_col)
|
||||
{
|
||||
if (line_length == (front_space+back_space+4))
|
||||
{
|
||||
for (int i = 0; i < front_space+4; i++) clog << " ";
|
||||
clog << segment << " ";
|
||||
line_length += (segment.length()+1);
|
||||
}
|
||||
else
|
||||
{
|
||||
clog << segment << " ";
|
||||
line_length += (segment.length()+1);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
clog << endl;
|
||||
for (int i = 0; i < front_space+9; i++) clog << " ";
|
||||
clog << segment << " ";
|
||||
line_length = (segment.length()+10+front_space+back_space);
|
||||
}
|
||||
}
|
||||
clog << endl;
|
||||
|
||||
if (!options[i].sec_message.empty())
|
||||
{
|
||||
for (int j = 0; j < options[i].sec_message.size(); j++)
|
||||
{
|
||||
ss_message.clear(); ss_message.str(options[i].sec_message[j]);
|
||||
|
||||
line_length = front_space + back_space + 9;
|
||||
while(ss_message >> segment)
|
||||
{
|
||||
if ((line_length+segment.length()+1) <= w.ws_col)
|
||||
{
|
||||
if (line_length == (front_space+back_space+9))
|
||||
{
|
||||
for (int i = 0; i < front_space+9; i++) clog << " ";
|
||||
clog << segment << " ";
|
||||
line_length += (segment.length()+1);
|
||||
}
|
||||
else
|
||||
{
|
||||
clog << segment << " ";
|
||||
line_length += (segment.length()+1);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
clog << endl;
|
||||
for (int i = 0; i < front_space+13; i++) clog << " ";
|
||||
clog << segment << " ";
|
||||
line_length = (segment.length()+14+front_space+back_space);
|
||||
}
|
||||
}
|
||||
clog << endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (!examples.empty())
|
||||
{
|
||||
for (int i = 0; i < front_space; i++) clog << " ";
|
||||
clog << "Examples:" << endl;
|
||||
for (int i = 0; i < examples.size(); i++)
|
||||
{
|
||||
ss_message.clear(); ss_message.str(examples[i]);
|
||||
|
||||
line_length = front_space + back_space + 4;
|
||||
while(ss_message >> segment)
|
||||
{
|
||||
if ((line_length+segment.length()+1) <= w.ws_col)
|
||||
{
|
||||
if (line_length == (front_space+back_space+4))
|
||||
{
|
||||
for (int i = 0; i < front_space+4; i++) clog << " ";
|
||||
clog << segment << " ";
|
||||
line_length += (segment.length()+1);
|
||||
}
|
||||
else
|
||||
{
|
||||
clog << segment << " ";
|
||||
line_length += (segment.length()+1);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
clog << endl;
|
||||
for (int i = 0; i < front_space+9; i++) clog << " ";
|
||||
clog << segment << " ";
|
||||
line_length = (segment.length()+10+front_space+back_space);
|
||||
}
|
||||
}
|
||||
clog << endl;
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
@ -1,56 +0,0 @@
|
||||
#ifndef _DISPHELP_H
|
||||
#define _DISPHELP_H
|
||||
#include <iostream>
|
||||
#include <sstream>
|
||||
#include <fstream>
|
||||
#include <string.h>
|
||||
#include <stdlib.h>
|
||||
#include <stdio.h>
|
||||
#include <unistd.h>
|
||||
#include <iomanip>
|
||||
#include <sys/ioctl.h>
|
||||
#include "vector"
|
||||
|
||||
using namespace std;
|
||||
|
||||
typedef vector<string> strArray;
|
||||
|
||||
struct option
|
||||
{
|
||||
string flag_s,flag_l;
|
||||
string message;
|
||||
strArray sec_message;
|
||||
option()
|
||||
{
|
||||
flag_s = flag_l = message = "";
|
||||
}
|
||||
};
|
||||
typedef vector<option> opArray;
|
||||
|
||||
class DispHelp
|
||||
{
|
||||
public:
|
||||
DispHelp(){
|
||||
front_space = 0;
|
||||
back_space = 10;
|
||||
ex_name = "Execuable";
|
||||
version = "0.0.1";
|
||||
descript = "Brief information about this command.";
|
||||
author = "Author's information.";
|
||||
}
|
||||
~DispHelp(){}
|
||||
void addHeadInfo(string,string,string,string);
|
||||
void addUsage(string);
|
||||
void addOption(string,string,string lflag = "");
|
||||
void addOptionSec(string,int index = -1);
|
||||
void addExample(string);
|
||||
void changeLayerOut(int,int);
|
||||
void show();
|
||||
private:
|
||||
string ex_name,version,descript,author;
|
||||
int front_space,back_space;
|
||||
opArray options;
|
||||
strArray examples;
|
||||
strArray usages;
|
||||
};
|
||||
#endif
|
||||
@ -1,6 +1,6 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::ForwardDeltaT(char* noise_level,char* mag_para){
|
||||
int GM3D::ForwardDeltaT(const char* noise_level,const char* mag_para){
|
||||
int i,j;
|
||||
double I0,A0,I,A;
|
||||
double k1,k2,k3,k4,k5,k6;
|
||||
@ -21,10 +21,10 @@ int GM3D::ForwardDeltaT(char* noise_level,char* mag_para){
|
||||
I0 = I = 90; A0 = A = 0;
|
||||
}
|
||||
else{
|
||||
I0=I0*Pi/180;
|
||||
A0=A0*Pi/180;
|
||||
I=I*Pi/180;
|
||||
A=A*Pi/180;
|
||||
I0=I0*GCTL_Pi/180;
|
||||
A0=A0*GCTL_Pi/180;
|
||||
I=I*GCTL_Pi/180;
|
||||
A=A*GCTL_Pi/180;
|
||||
|
||||
k1=cos(I0)*sin(A0)*sin(I)+sin(I0)*cos(I)*sin(A);
|
||||
k2=cos(I0)*cos(A0)*sin(I)+sin(I0)*cos(I)*cos(A);
|
||||
@ -38,18 +38,18 @@ int GM3D::ForwardDeltaT(char* noise_level,char* mag_para){
|
||||
default_random_engine generator;
|
||||
normal_distribution<double> dist(noise_mean, noise_dev);
|
||||
|
||||
_1dArray temp_obs(model_num_);
|
||||
gctl::_1d_vector temp_obs(model_num_);
|
||||
|
||||
ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
|
||||
gctl::progress_bar bar(obs_num_,"Forward modeling");
|
||||
for (i = 0; i < obs_num_; i++){
|
||||
bar->Progressed(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 < model_num_; j++){
|
||||
if (fabs(forward_model_[j]) > ZERO){
|
||||
if (fabs(forward_model_[j]) > GCTL_ZERO){
|
||||
x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
|
||||
y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
|
||||
z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
|
||||
@ -113,7 +113,7 @@ int GM3D::ForwardDeltaT(char* noise_level,char* mag_para){
|
||||
+k4*atan(((x1-obs_p_[i].x)*(y1-obs_p_[i].y))/w111)+k5*atan(((x1-obs_p_[i].x)*(y1-obs_p_[i].y))/v111)
|
||||
+k6*atan(((x1-obs_p_[i].x)*(y1-obs_p_[i].y))/(R111*(z1-obs_p_[i].z)));
|
||||
|
||||
temp_obs[j] = T0/(4*Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
temp_obs[j] = GCTL_T0/(4*GCTL_Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::ForwardDeltaTx(char* noise_level,char* mag_para){
|
||||
int GM3D::ForwardDeltaTx(const char* noise_level,const char* mag_para){
|
||||
int i,j;
|
||||
double I0,A0,I,A;
|
||||
double k1,k2,k3,k4,k5,k6;
|
||||
@ -21,10 +21,10 @@ int GM3D::ForwardDeltaTx(char* noise_level,char* mag_para){
|
||||
I0 = I = 90; A0 = A = 0;
|
||||
}
|
||||
else{
|
||||
I0=I0*Pi/180;
|
||||
A0=A0*Pi/180;
|
||||
I=I*Pi/180;
|
||||
A=A*Pi/180;
|
||||
I0=I0*GCTL_Pi/180;
|
||||
A0=A0*GCTL_Pi/180;
|
||||
I=I*GCTL_Pi/180;
|
||||
A=A*GCTL_Pi/180;
|
||||
|
||||
k1=cos(I0)*sin(A0)*sin(I)+sin(I0)*cos(I)*sin(A);
|
||||
k2=cos(I0)*cos(A0)*sin(I)+sin(I0)*cos(I)*cos(A);
|
||||
@ -38,18 +38,18 @@ int GM3D::ForwardDeltaTx(char* noise_level,char* mag_para){
|
||||
default_random_engine generator;
|
||||
normal_distribution<double> dist(noise_mean, noise_dev);
|
||||
|
||||
_1dArray temp_obs(model_num_);
|
||||
gctl::_1d_vector temp_obs(model_num_);
|
||||
|
||||
ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
|
||||
gctl::progress_bar bar(obs_num_,"Forward modeling");
|
||||
for (i = 0; i < obs_num_; i++){
|
||||
bar->Progressed(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 < model_num_; j++){
|
||||
if (fabs(forward_model_[j]) > ZERO){
|
||||
if (fabs(forward_model_[j]) > GCTL_ZERO){
|
||||
x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
|
||||
y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
|
||||
z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
|
||||
@ -90,7 +90,7 @@ int GM3D::ForwardDeltaTx(char* noise_level,char* mag_para){
|
||||
G211=-k1*pow(R211,-1)-k2*(x2-obs_p_[i].x)/(R211*(R211+y1-obs_p_[i].y))-k3*(x2-obs_p_[i].x)/(R211*(R211+z1-obs_p_[i].z))+k4*w211-k5*(y1-obs_p_[i].y)/(R211*(R211+z1-obs_p_[i].z))-k6*v211;
|
||||
G111=-k1*pow(R111,-1)-k2*(x1-obs_p_[i].x)/(R111*(R111+y1-obs_p_[i].y))-k3*(x1-obs_p_[i].x)/(R111*(R111+z1-obs_p_[i].z))+k4*w111-k5*(y1-obs_p_[i].y)/(R111*(R111+z1-obs_p_[i].z))-k6*v111;
|
||||
|
||||
temp_obs[j] = T0/(4*Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
temp_obs[j] = GCTL_T0/(4*GCTL_Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::ForwardDeltaTy(char* noise_level,char* mag_para){
|
||||
int GM3D::ForwardDeltaTy(const char* noise_level,const char* mag_para){
|
||||
int i,j;
|
||||
double I0,A0,I,A;
|
||||
double k1,k2,k3,k4,k5,k6;
|
||||
@ -21,10 +21,10 @@ int GM3D::ForwardDeltaTy(char* noise_level,char* mag_para){
|
||||
I0 = I = 90; A0 = A = 0;
|
||||
}
|
||||
else{
|
||||
I0=I0*Pi/180;
|
||||
A0=A0*Pi/180;
|
||||
I=I*Pi/180;
|
||||
A=A*Pi/180;
|
||||
I0=I0*GCTL_Pi/180;
|
||||
A0=A0*GCTL_Pi/180;
|
||||
I=I*GCTL_Pi/180;
|
||||
A=A*GCTL_Pi/180;
|
||||
|
||||
k1=cos(I0)*sin(A0)*sin(I)+sin(I0)*cos(I)*sin(A);
|
||||
k2=cos(I0)*cos(A0)*sin(I)+sin(I0)*cos(I)*cos(A);
|
||||
@ -38,18 +38,18 @@ int GM3D::ForwardDeltaTy(char* noise_level,char* mag_para){
|
||||
default_random_engine generator;
|
||||
normal_distribution<double> dist(noise_mean, noise_dev);
|
||||
|
||||
_1dArray temp_obs(model_num_);
|
||||
gctl::_1d_vector temp_obs(model_num_);
|
||||
|
||||
ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
|
||||
gctl::progress_bar bar(obs_num_,"Forward modeling");
|
||||
for (i = 0; i < obs_num_; i++){
|
||||
bar->Progressed(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 < model_num_; j++){
|
||||
if (fabs(forward_model_[j]) > ZERO){
|
||||
if (fabs(forward_model_[j]) > GCTL_ZERO){
|
||||
x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
|
||||
y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
|
||||
z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
|
||||
@ -90,7 +90,7 @@ int GM3D::ForwardDeltaTy(char* noise_level,char* mag_para){
|
||||
G211=-k1*(y1-obs_p_[i].y)/(R211*(R211+x2-obs_p_[i].x))-k2*pow(R211,-1)-k3*(y1-obs_p_[i].y)/(R211*(R211+z1-obs_p_[i].z))-k4*(x2-obs_p_[i].x)/(R211*(R211+z1-obs_p_[i].z))+k5*w211-k6*v211;
|
||||
G111=-k1*(y1-obs_p_[i].y)/(R111*(R111+x1-obs_p_[i].x))-k2*pow(R111,-1)-k3*(y1-obs_p_[i].y)/(R111*(R111+z1-obs_p_[i].z))-k4*(x1-obs_p_[i].x)/(R111*(R111+z1-obs_p_[i].z))+k5*w111-k6*v111;
|
||||
|
||||
temp_obs[j] = T0/(4*Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
temp_obs[j] = GCTL_T0/(4*GCTL_Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::ForwardDeltaTz(char* noise_level,char* mag_para){
|
||||
int GM3D::ForwardDeltaTz(const char* noise_level,const char* mag_para){
|
||||
int i,j;
|
||||
double I0,A0,I,A;
|
||||
double k1,k2,k3,k4,k5,k6;
|
||||
@ -22,10 +22,10 @@ int GM3D::ForwardDeltaTz(char* noise_level,char* mag_para){
|
||||
I0 = I = 90; A0 = A = 0;
|
||||
}
|
||||
else{
|
||||
I0=I0*Pi/180;
|
||||
A0=A0*Pi/180;
|
||||
I=I*Pi/180;
|
||||
A=A*Pi/180;
|
||||
I0=I0*GCTL_Pi/180;
|
||||
A0=A0*GCTL_Pi/180;
|
||||
I=I*GCTL_Pi/180;
|
||||
A=A*GCTL_Pi/180;
|
||||
|
||||
k1=cos(I0)*sin(A0)*sin(I)+sin(I0)*cos(I)*sin(A);
|
||||
k2=cos(I0)*cos(A0)*sin(I)+sin(I0)*cos(I)*cos(A);
|
||||
@ -39,11 +39,11 @@ int GM3D::ForwardDeltaTz(char* noise_level,char* mag_para){
|
||||
default_random_engine generator;
|
||||
normal_distribution<double> dist(noise_mean, noise_dev);
|
||||
|
||||
_1dArray temp_obs(model_num_);
|
||||
gctl::_1d_vector temp_obs(model_num_);
|
||||
|
||||
ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
|
||||
gctl::progress_bar bar(obs_num_,"Forward modeling");
|
||||
for (i = 0; i < obs_num_; i++){
|
||||
bar->Progressed(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,\
|
||||
t111,t112,t121,t122,t211,t212,t221,t222,\
|
||||
@ -51,7 +51,7 @@ int GM3D::ForwardDeltaTz(char* noise_level,char* mag_para){
|
||||
v111,v112,v121,v122,v211,v212,v221,v222,\
|
||||
x1,x2,y1,y2,z1,z2) shared(i) schedule(guided)
|
||||
for (j = 0; j < model_num_; j++){
|
||||
if (fabs(forward_model_[j]) > ZERO){
|
||||
if (fabs(forward_model_[j]) > GCTL_ZERO){
|
||||
x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
|
||||
y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
|
||||
z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
|
||||
@ -101,7 +101,7 @@ int GM3D::ForwardDeltaTz(char* noise_level,char* mag_para){
|
||||
G211=-k1*(z1-obs_p_[i].z)/(R211*(R211+x2-obs_p_[i].x))-k2*(z1-obs_p_[i].z)/(R211*(R211+y1-obs_p_[i].y))-k3/R211+k4*t211+k5*v211+k6*w211;
|
||||
G111=-k1*(z1-obs_p_[i].z)/(R111*(R111+x1-obs_p_[i].x))-k2*(z1-obs_p_[i].z)/(R111*(R111+y1-obs_p_[i].y))-k3/R111+k4*t111+k5*v111+k6*w111;
|
||||
|
||||
temp_obs[j] = T0/(4*Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
temp_obs[j] = GCTL_T0/(4*GCTL_Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::ForwardHax(char* noise_level,char* mag_para){
|
||||
int GM3D::ForwardHax(const char* noise_level,const char* mag_para){
|
||||
int i,j;
|
||||
double I0,A0,I,A;
|
||||
double Alpha,Beta,Gamma;
|
||||
@ -21,8 +21,8 @@ int GM3D::ForwardHax(char* noise_level,char* mag_para){
|
||||
I = 90; A = 0;
|
||||
}
|
||||
else{
|
||||
I=I*Pi/180;
|
||||
A=A*Pi/180;
|
||||
I=I*GCTL_Pi/180;
|
||||
A=A*GCTL_Pi/180;
|
||||
|
||||
Alpha=cos(I)*cos(A);
|
||||
Beta=cos(I)*sin(A);
|
||||
@ -33,18 +33,18 @@ int GM3D::ForwardHax(char* noise_level,char* mag_para){
|
||||
default_random_engine generator;
|
||||
normal_distribution<double> dist(noise_mean, noise_dev);
|
||||
|
||||
_1dArray temp_obs(model_num_);
|
||||
gctl::_1d_vector temp_obs(model_num_);
|
||||
|
||||
ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
|
||||
gctl::progress_bar bar(obs_num_,"Forward modeling");
|
||||
for (i = 0; i < obs_num_; i++){
|
||||
bar->Progressed(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_1,w112_1,w121_1,w122_1,w211_1,w212_1,w221_1,w222_1,\
|
||||
w111_2,w112_2,w121_2,w122_2,w211_2,w212_2,w221_2,w222_2,\
|
||||
x1,x2,y1,y2,z1,z2) shared(i) schedule(guided)
|
||||
for (j = 0; j < model_num_; j++){
|
||||
if (fabs(forward_model_[j]) > ZERO){
|
||||
if (fabs(forward_model_[j]) > GCTL_ZERO){
|
||||
x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
|
||||
y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
|
||||
z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
|
||||
@ -85,7 +85,7 @@ int GM3D::ForwardHax(char* noise_level,char* mag_para){
|
||||
G211=Alpha*atan2(w211_1,w211_2)+Beta*log(R211+z1-obs_p_[i].z)+Gamma*log(R211+y1-obs_p_[i].y);
|
||||
G111=Alpha*atan2(w111_1,w111_2)+Beta*log(R111+z1-obs_p_[i].z)+Gamma*log(R111+y1-obs_p_[i].y);
|
||||
|
||||
temp_obs[j] = T0/(4*Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
temp_obs[j] = GCTL_T0/(4*GCTL_Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::ForwardHay(char* noise_level,char* mag_para){
|
||||
int GM3D::ForwardHay(const char* noise_level,const char* mag_para){
|
||||
int i,j;
|
||||
double I0,A0,I,A;
|
||||
double Alpha,Beta,Gamma;
|
||||
@ -21,8 +21,8 @@ int GM3D::ForwardHay(char* noise_level,char* mag_para){
|
||||
I = 90; A = 0;
|
||||
}
|
||||
else{
|
||||
I=I*Pi/180;
|
||||
A=A*Pi/180;
|
||||
I=I*GCTL_Pi/180;
|
||||
A=A*GCTL_Pi/180;
|
||||
|
||||
Alpha=cos(I)*cos(A);
|
||||
Beta=cos(I)*sin(A);
|
||||
@ -33,18 +33,18 @@ int GM3D::ForwardHay(char* noise_level,char* mag_para){
|
||||
default_random_engine generator;
|
||||
normal_distribution<double> dist(noise_mean, noise_dev);
|
||||
|
||||
_1dArray temp_obs(model_num_);
|
||||
gctl::_1d_vector temp_obs(model_num_);
|
||||
|
||||
ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
|
||||
gctl::progress_bar bar(obs_num_,"Forward modeling");
|
||||
for (i = 0; i < obs_num_; i++){
|
||||
bar->Progressed(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_1,w112_1,w121_1,w122_1,w211_1,w212_1,w221_1,w222_1,\
|
||||
w111_2,w112_2,w121_2,w122_2,w211_2,w212_2,w221_2,w222_2,\
|
||||
x1,x2,y1,y2,z1,z2) shared(i) schedule(guided)
|
||||
for (j = 0; j < model_num_; j++){
|
||||
if (fabs(forward_model_[j]) > ZERO){
|
||||
if (fabs(forward_model_[j]) > GCTL_ZERO){
|
||||
x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
|
||||
y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
|
||||
z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
|
||||
@ -85,7 +85,7 @@ int GM3D::ForwardHay(char* noise_level,char* mag_para){
|
||||
G211=Beta*atan2(w211_1,w211_2)+Alpha*log(R211+z1-obs_p_[i].z)+Gamma*log(R211+x2-obs_p_[i].x);
|
||||
G111=Beta*atan2(w111_1,w111_2)+Alpha*log(R111+z1-obs_p_[i].z)+Gamma*log(R111+x1-obs_p_[i].x);
|
||||
|
||||
temp_obs[j] = T0/(4*Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
temp_obs[j] = GCTL_T0/(4*GCTL_Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::ForwardVz(char* noise_level){
|
||||
int GM3D::ForwardVz(const char* noise_level){
|
||||
int i,j;
|
||||
double x1,x2,y1,y2,z1,z2;
|
||||
double R222,R122,R212,R112,R221,R121,R211,R111;
|
||||
@ -16,14 +16,14 @@ int GM3D::ForwardVz(char* noise_level){
|
||||
default_random_engine generator;
|
||||
normal_distribution<double> dist(noise_mean, noise_dev);
|
||||
|
||||
_1dArray temp_obs(model_num_);
|
||||
gctl::_1d_vector temp_obs(model_num_);
|
||||
|
||||
ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
|
||||
gctl::progress_bar bar(obs_num_,"Forward modeling");
|
||||
for (i = 0; i < obs_num_; i++){
|
||||
bar->Progressed(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,x1,x2,y1,y2,z1,z2) shared(i) schedule(guided)
|
||||
for (j = 0; j < model_num_; j++){
|
||||
if (fabs(forward_model_[j]) > ZERO){
|
||||
if (fabs(forward_model_[j]) > GCTL_ZERO){
|
||||
x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
|
||||
y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
|
||||
z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
|
||||
@ -37,16 +37,16 @@ int GM3D::ForwardVz(char* noise_level){
|
||||
R211=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
|
||||
R111=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
|
||||
|
||||
G222=(x2-obs_p_[i].x)*log((y2-obs_p_[i].y)+R222)+(y2-obs_p_[i].y)*log((x2-obs_p_[i].x)+R222)+(z2-obs_p_[i].z)*arctg((z2-obs_p_[i].z)*R222/(x2-obs_p_[i].x)/(y2-obs_p_[i].y));
|
||||
G122=(x1-obs_p_[i].x)*log((y2-obs_p_[i].y)+R122)+(y2-obs_p_[i].y)*log((x1-obs_p_[i].x)+R122)+(z2-obs_p_[i].z)*arctg((z2-obs_p_[i].z)*R122/(x1-obs_p_[i].x)/(y2-obs_p_[i].y));
|
||||
G212=(x2-obs_p_[i].x)*log((y1-obs_p_[i].y)+R212)+(y1-obs_p_[i].y)*log((x2-obs_p_[i].x)+R212)+(z2-obs_p_[i].z)*arctg((z2-obs_p_[i].z)*R212/(x2-obs_p_[i].x)/(y1-obs_p_[i].y));
|
||||
G112=(x1-obs_p_[i].x)*log((y1-obs_p_[i].y)+R112)+(y1-obs_p_[i].y)*log((x1-obs_p_[i].x)+R112)+(z2-obs_p_[i].z)*arctg((z2-obs_p_[i].z)*R112/(x1-obs_p_[i].x)/(y1-obs_p_[i].y));
|
||||
G221=(x2-obs_p_[i].x)*log((y2-obs_p_[i].y)+R221)+(y2-obs_p_[i].y)*log((x2-obs_p_[i].x)+R221)+(z1-obs_p_[i].z)*arctg((z1-obs_p_[i].z)*R221/(x2-obs_p_[i].x)/(y2-obs_p_[i].y));
|
||||
G121=(x1-obs_p_[i].x)*log((y2-obs_p_[i].y)+R121)+(y2-obs_p_[i].y)*log((x1-obs_p_[i].x)+R121)+(z1-obs_p_[i].z)*arctg((z1-obs_p_[i].z)*R121/(x1-obs_p_[i].x)/(y2-obs_p_[i].y));
|
||||
G211=(x2-obs_p_[i].x)*log((y1-obs_p_[i].y)+R211)+(y1-obs_p_[i].y)*log((x2-obs_p_[i].x)+R211)+(z1-obs_p_[i].z)*arctg((z1-obs_p_[i].z)*R211/(x2-obs_p_[i].x)/(y1-obs_p_[i].y));
|
||||
G111=(x1-obs_p_[i].x)*log((y1-obs_p_[i].y)+R111)+(y1-obs_p_[i].y)*log((x1-obs_p_[i].x)+R111)+(z1-obs_p_[i].z)*arctg((z1-obs_p_[i].z)*R111/(x1-obs_p_[i].x)/(y1-obs_p_[i].y));
|
||||
G222=(x2-obs_p_[i].x)*log((y2-obs_p_[i].y)+R222)+(y2-obs_p_[i].y)*log((x2-obs_p_[i].x)+R222)+(z2-obs_p_[i].z)*gctl::arctg((z2-obs_p_[i].z)*R222/(x2-obs_p_[i].x)/(y2-obs_p_[i].y));
|
||||
G122=(x1-obs_p_[i].x)*log((y2-obs_p_[i].y)+R122)+(y2-obs_p_[i].y)*log((x1-obs_p_[i].x)+R122)+(z2-obs_p_[i].z)*gctl::arctg((z2-obs_p_[i].z)*R122/(x1-obs_p_[i].x)/(y2-obs_p_[i].y));
|
||||
G212=(x2-obs_p_[i].x)*log((y1-obs_p_[i].y)+R212)+(y1-obs_p_[i].y)*log((x2-obs_p_[i].x)+R212)+(z2-obs_p_[i].z)*gctl::arctg((z2-obs_p_[i].z)*R212/(x2-obs_p_[i].x)/(y1-obs_p_[i].y));
|
||||
G112=(x1-obs_p_[i].x)*log((y1-obs_p_[i].y)+R112)+(y1-obs_p_[i].y)*log((x1-obs_p_[i].x)+R112)+(z2-obs_p_[i].z)*gctl::arctg((z2-obs_p_[i].z)*R112/(x1-obs_p_[i].x)/(y1-obs_p_[i].y));
|
||||
G221=(x2-obs_p_[i].x)*log((y2-obs_p_[i].y)+R221)+(y2-obs_p_[i].y)*log((x2-obs_p_[i].x)+R221)+(z1-obs_p_[i].z)*gctl::arctg((z1-obs_p_[i].z)*R221/(x2-obs_p_[i].x)/(y2-obs_p_[i].y));
|
||||
G121=(x1-obs_p_[i].x)*log((y2-obs_p_[i].y)+R121)+(y2-obs_p_[i].y)*log((x1-obs_p_[i].x)+R121)+(z1-obs_p_[i].z)*gctl::arctg((z1-obs_p_[i].z)*R121/(x1-obs_p_[i].x)/(y2-obs_p_[i].y));
|
||||
G211=(x2-obs_p_[i].x)*log((y1-obs_p_[i].y)+R211)+(y1-obs_p_[i].y)*log((x2-obs_p_[i].x)+R211)+(z1-obs_p_[i].z)*gctl::arctg((z1-obs_p_[i].z)*R211/(x2-obs_p_[i].x)/(y1-obs_p_[i].y));
|
||||
G111=(x1-obs_p_[i].x)*log((y1-obs_p_[i].y)+R111)+(y1-obs_p_[i].y)*log((x1-obs_p_[i].x)+R111)+(z1-obs_p_[i].z)*gctl::arctg((z1-obs_p_[i].z)*R111/(x1-obs_p_[i].x)/(y1-obs_p_[i].y));
|
||||
|
||||
temp_obs[j] = -1.0*G0*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
temp_obs[j] = -1.0e+8*GCTL_G0*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::ForwardVzx(char* noise_level){
|
||||
int GM3D::ForwardVzx(const char* noise_level){
|
||||
int i,j;
|
||||
double x1,x2,y1,y2,z1,z2;
|
||||
double R222,R122,R212,R112,R221,R121,R211,R111;
|
||||
@ -15,14 +15,14 @@ int GM3D::ForwardVzx(char* noise_level){
|
||||
default_random_engine generator;
|
||||
normal_distribution<double> dist(noise_mean, noise_dev);
|
||||
|
||||
_1dArray temp_obs(model_num_);
|
||||
gctl::_1d_vector temp_obs(model_num_);
|
||||
|
||||
ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
|
||||
gctl::progress_bar bar(obs_num_,"Forward modeling");
|
||||
for (i = 0; i < obs_num_; i++){
|
||||
bar->Progressed(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,x1,x2,y1,y2,z1,z2) shared(i) schedule(guided)
|
||||
for (j = 0; j < model_num_; j++){
|
||||
if (fabs(forward_model_[j]) > ZERO){
|
||||
if (fabs(forward_model_[j]) > GCTL_ZERO){
|
||||
x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
|
||||
y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
|
||||
z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
|
||||
@ -45,7 +45,7 @@ int GM3D::ForwardVzx(char* noise_level){
|
||||
G211=log((y1-obs_p_[i].y)+R211);
|
||||
G111=log((y1-obs_p_[i].y)+R111);
|
||||
|
||||
temp_obs[j] = G0*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
temp_obs[j] = 1.0e+8*GCTL_G0*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::ForwardVzy(char* noise_level){
|
||||
int GM3D::ForwardVzy(const char* noise_level){
|
||||
int i,j;
|
||||
double x1,x2,y1,y2,z1,z2;
|
||||
double R222,R122,R212,R112,R221,R121,R211,R111;
|
||||
@ -15,14 +15,14 @@ int GM3D::ForwardVzy(char* noise_level){
|
||||
default_random_engine generator;
|
||||
normal_distribution<double> dist(noise_mean, noise_dev);
|
||||
|
||||
_1dArray temp_obs(model_num_);
|
||||
gctl::_1d_vector temp_obs(model_num_);
|
||||
|
||||
ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
|
||||
gctl::progress_bar bar(obs_num_,"Forward modeling");
|
||||
for (i = 0; i < obs_num_; i++){
|
||||
bar->Progressed(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,x1,x2,y1,y2,z1,z2) shared(i) schedule(guided)
|
||||
for (j = 0; j < model_num_; j++){
|
||||
if (fabs(forward_model_[j]) > ZERO){
|
||||
if (fabs(forward_model_[j]) > GCTL_ZERO){
|
||||
x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
|
||||
y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
|
||||
z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
|
||||
@ -45,7 +45,7 @@ int GM3D::ForwardVzy(char* noise_level){
|
||||
G211=log((x2-obs_p_[i].x)+R211);
|
||||
G111=log((x1-obs_p_[i].x)+R111);
|
||||
|
||||
temp_obs[j] = G0*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
temp_obs[j] = 1.0e+8*GCTL_G0*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::ForwardVzz(char* noise_level){
|
||||
int GM3D::ForwardVzz(const char* noise_level){
|
||||
int i,j;
|
||||
double x1,x2,y1,y2,z1,z2;
|
||||
double R222,R122,R212,R112,R221,R121,R211,R111;
|
||||
@ -15,14 +15,14 @@ int GM3D::ForwardVzz(char* noise_level){
|
||||
default_random_engine generator;
|
||||
normal_distribution<double> dist(noise_mean, noise_dev);
|
||||
|
||||
_1dArray temp_obs(model_num_);
|
||||
gctl::_1d_vector temp_obs(model_num_);
|
||||
|
||||
ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
|
||||
gctl::progress_bar bar(obs_num_,"Forward modeling");
|
||||
for (i = 0; i < obs_num_; i++){
|
||||
bar->Progressed(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,x1,x2,y1,y2,z1,z2) shared(i) schedule(guided)
|
||||
for (j = 0; j < model_num_; j++){
|
||||
if (fabs(forward_model_[j]) > ZERO){
|
||||
if (fabs(forward_model_[j]) > GCTL_ZERO){
|
||||
x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
|
||||
y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
|
||||
z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
|
||||
@ -54,7 +54,7 @@ int GM3D::ForwardVzz(char* noise_level){
|
||||
G211=atan((x2-obs_p_[i].x)*(y1-obs_p_[i].y)/(R211*(z1-obs_p_[i].z)));
|
||||
G111=atan((x1-obs_p_[i].x)*(y1-obs_p_[i].y)/(R111*(z1-obs_p_[i].z)));
|
||||
|
||||
temp_obs[j] = -1.0*G0*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
temp_obs[j] = -1.0e+8*GCTL_G0*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::ForwardZa(char* noise_level,char* mag_para){
|
||||
int GM3D::ForwardZa(const char* noise_level,const char* mag_para){
|
||||
int i,j;
|
||||
double I0,A0,I,A;
|
||||
double Alpha,Beta,Gamma;
|
||||
@ -21,8 +21,8 @@ int GM3D::ForwardZa(char* noise_level,char* mag_para){
|
||||
I = 90; A = 0;
|
||||
}
|
||||
else{
|
||||
I=I*Pi/180;
|
||||
A=A*Pi/180;
|
||||
I=I*GCTL_Pi/180;
|
||||
A=A*GCTL_Pi/180;
|
||||
|
||||
Alpha=cos(I)*cos(A);
|
||||
Beta=cos(I)*sin(A);
|
||||
@ -33,18 +33,18 @@ int GM3D::ForwardZa(char* noise_level,char* mag_para){
|
||||
default_random_engine generator;
|
||||
normal_distribution<double> dist(noise_mean, noise_dev);
|
||||
|
||||
_1dArray temp_obs(model_num_);
|
||||
gctl::_1d_vector temp_obs(model_num_);
|
||||
|
||||
ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
|
||||
gctl::progress_bar bar(obs_num_,"Forward modeling");
|
||||
for (i = 0; i < obs_num_; i++){
|
||||
bar->Progressed(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_1,w112_1,w121_1,w122_1,w211_1,w212_1,w221_1,w222_1,\
|
||||
w111_2,w112_2,w121_2,w122_2,w211_2,w212_2,w221_2,w222_2,\
|
||||
x1,x2,y1,y2,z1,z2) shared(i) schedule(guided)
|
||||
for (j = 0; j < model_num_; j++){
|
||||
if (fabs(forward_model_[j]) > ZERO){
|
||||
if (fabs(forward_model_[j]) > GCTL_ZERO){
|
||||
x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
|
||||
y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
|
||||
z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
|
||||
@ -85,7 +85,7 @@ int GM3D::ForwardZa(char* noise_level,char* mag_para){
|
||||
G211=-Gamma*atan2(w211_1,w211_2)+Alpha*log(R211+y1-obs_p_[i].y)+Beta*log(R211+x2-obs_p_[i].x);
|
||||
G111=-Gamma*atan2(w111_1,w111_2)+Alpha*log(R111+y1-obs_p_[i].y)+Beta*log(R111+x1-obs_p_[i].x);
|
||||
|
||||
temp_obs[j] = T0/(4*Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
temp_obs[j] = GCTL_T0/(4*GCTL_Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
139
src/gm3d.h
139
src/gm3d.h
@ -1,58 +1,119 @@
|
||||
#ifndef _GM3D_H
|
||||
#define _GM3D_H
|
||||
#include "head_func.h"
|
||||
#include "progress_bar.h"
|
||||
|
||||
#include "stdio.h"
|
||||
#include "stdlib.h"
|
||||
#include "unistd.h"
|
||||
#include "map"
|
||||
#include "algorithm"
|
||||
#include "ctime"
|
||||
#include "omp.h"
|
||||
#include "random"
|
||||
|
||||
// add gctl head files here
|
||||
#include "gctl/gctl_core.h"
|
||||
#include "gctl/gctl_utility.h"
|
||||
|
||||
using namespace std;
|
||||
|
||||
// define error code here
|
||||
#define GM3D_INVALID_MESH_PARAMETER "Invalid mesh parameter."
|
||||
#define GM3D_MODEL_INIT_FAIL "Fail to read models' parameter."
|
||||
#define GM3D_FORWARD_MODEL_INITIALIZED "Forward model initialized from file."
|
||||
#define GM3D_INVALID_OBS_PARAMETER "Invalid observation parameter."
|
||||
#define GM3D_NOMODEL_ADDED "No model changed."
|
||||
#define GM3D_INVALID_MODEL_TYPE "Wrong parameter for model's type."
|
||||
#define GM3D_INVALID_MODEL_PARAMETER "Invalid model parameter."
|
||||
#define GM3D_INVALID_MODEL_VAL_TYPE "Wrong parameter for model's value type."
|
||||
|
||||
//数据结构
|
||||
typedef map<int,int> int2intMap; ///< 整型到整形的映射
|
||||
|
||||
//模型块体参数
|
||||
struct modelist
|
||||
{
|
||||
std::string mod_type;
|
||||
std::string val_type;
|
||||
std::string mod_para;
|
||||
double mod_value;
|
||||
};
|
||||
typedef vector<modelist> modelistVector;
|
||||
//直角坐标系顶点 比坐标点位只多一个序号
|
||||
struct vertex : public gctl::point3d_c
|
||||
{
|
||||
int id = -1;
|
||||
};
|
||||
typedef vector<vertex> vertexVector;
|
||||
typedef map<string, vertex> str2vertMap;
|
||||
//观测点
|
||||
struct obspoint : public vertex
|
||||
{
|
||||
double val = GCTL_BDL_MAX; double dev = GCTL_BDL_MAX;
|
||||
};
|
||||
typedef vector<obspoint> obspointVector;
|
||||
//块体
|
||||
struct cube
|
||||
{
|
||||
vertex cen;
|
||||
int ids[8] = {-1,-1,-1,-1,-1,-1,-1,-1};
|
||||
double dx = GCTL_BDL_MAX; double dy = GCTL_BDL_MAX; double dz = GCTL_BDL_MAX;
|
||||
};
|
||||
typedef vector<cube> cubeVector;
|
||||
|
||||
class GM3D{
|
||||
public:
|
||||
GM3D(){}
|
||||
~GM3D(){}
|
||||
int BuildRegularGrid(char*); //初始化反演模型空间
|
||||
int AddModels(char*); //读取模型块体参数文件
|
||||
int AddRegularBlock(modelist); //添加普通模型块体
|
||||
int AddTiltedBlock(modelist); //添加倾斜模型块体
|
||||
int AddSphereBlock(modelist); //添加球体椭球体块体
|
||||
int AddInterfaceBlock(modelist); //添加密度界面
|
||||
void BuildRegularGrid(std::string modspace_para); //初始化反演模型空间
|
||||
void AddModels(std::string modentity_file); //读取模型块体参数文件
|
||||
void AddRegularBlock(modelist); //添加普通模型块体
|
||||
void AddTiltedBlock(modelist); //添加倾斜模型块体
|
||||
void AddSphereBlock(modelist); //添加球体椭球体块体
|
||||
void AddInterfaceBlock(modelist); //添加密度界面
|
||||
//模型操作
|
||||
int ReadModel(char*,char*);
|
||||
int ReadModel(const char*,const char*);
|
||||
//输出模型
|
||||
int RegisteredOuput(bool); //注册输出的块体模型
|
||||
int OutMshFile(char*,string); //输出模型文件
|
||||
int OutNeighborFile(char*,char*); //输出模型块体或顶点的相邻关系 暂缓
|
||||
int OutMshFile(const char*,string); //输出模型文件
|
||||
int OutNeighborFile(const char*,const char*); //输出模型块体或顶点的相邻关系 暂缓
|
||||
//观测数据
|
||||
int InitObs(char*);
|
||||
int OutObs(char*);
|
||||
int InitObs(const char*);
|
||||
int OutObs(const char*);
|
||||
//正演函数
|
||||
int ForwardVz(char*);
|
||||
int ForwardVzx(char*);
|
||||
int ForwardVzy(char*);
|
||||
int ForwardVzz(char*);
|
||||
int ForwardDeltaT(char*,char*);
|
||||
int ForwardDeltaTx(char*,char*);
|
||||
int ForwardDeltaTy(char*,char*);
|
||||
int ForwardDeltaTz(char*,char*);
|
||||
int ForwardHax(char*,char*);
|
||||
int ForwardHay(char*,char*);
|
||||
int ForwardZa(char*,char*);
|
||||
int ForwardVz(const char*);
|
||||
int ForwardVzx(const char*);
|
||||
int ForwardVzy(const char*);
|
||||
int ForwardVzz(const char*);
|
||||
int ForwardDeltaT(const char*,const char*);
|
||||
int ForwardDeltaTx(const char*,const char*);
|
||||
int ForwardDeltaTy(const char*,const char*);
|
||||
int ForwardDeltaTz(const char*,const char*);
|
||||
int ForwardHax(const char*,const char*);
|
||||
int ForwardHay(const char*,const char*);
|
||||
int ForwardZa(const char*,const char*);
|
||||
//返回一个cpoint的位置id
|
||||
string point_id(vertex);
|
||||
public:
|
||||
char model_id_[1024];
|
||||
private:
|
||||
int obs_num_, model_num_, vert_num_;
|
||||
//正演数组
|
||||
obspointArray obs_p_;
|
||||
_2dArray input_models_;
|
||||
_1sArray input_model_names_;
|
||||
_1dArray forward_model_;
|
||||
obspointVector obs_p_;
|
||||
gctl::_2d_vector input_models_;
|
||||
gctl::_1s_vector input_model_names_;
|
||||
gctl::_1d_vector forward_model_;
|
||||
//模型数据
|
||||
cubeArray model_cube_;
|
||||
cpointArray model_vert_;
|
||||
_1dArray model_block_val_;
|
||||
modelistArray model_list_;
|
||||
_1iArray out_ele_ids_;
|
||||
_1iArray out_ele_data_ids_;
|
||||
_1iArray out_vert_ids_;
|
||||
_int2intMap vert_out_map_;
|
||||
_int2intMap ele_data_out_map_;
|
||||
cubeVector model_cube_;
|
||||
vertexVector model_vert_;
|
||||
gctl::_1d_vector model_block_val_;
|
||||
modelistVector model_list_;
|
||||
gctl::_1i_vector out_ele_ids_;
|
||||
gctl::_1i_vector out_ele_data_ids_;
|
||||
gctl::_1i_vector out_vert_ids_;
|
||||
int2intMap vert_out_map_;
|
||||
int2intMap ele_data_out_map_;
|
||||
|
||||
_2iArray model_vert_neighbor_; //暂缓
|
||||
_2iArray model_cube_neighbor_; //暂缓
|
||||
gctl::_2i_vector model_vert_neighbor_; //暂缓
|
||||
gctl::_2i_vector model_cube_neighbor_; //暂缓
|
||||
};
|
||||
#endif
|
||||
@ -1,95 +0,0 @@
|
||||
#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模长
|
||||
double modCpoint(cpoint v){
|
||||
return sqrt(v.x*v.x+v.y*v.y+v.z*v.z);
|
||||
}
|
||||
/*************************全局函数********************************/
|
||||
//正负分离的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
|
||||
stringstream str2ss(string s){
|
||||
stringstream sstr;
|
||||
sstr.str(""); sstr.clear(); sstr.str(s);
|
||||
return sstr;
|
||||
}
|
||||
//返回一个cpoint的位置id字符串
|
||||
string cpoint_id(cpoint c){
|
||||
string vert_id, mid_id;
|
||||
stringstream sstemp;
|
||||
sstemp.str(""); sstemp.clear(); sstemp<<setprecision(PRECISION)<<c.x;
|
||||
sstemp>>vert_id;
|
||||
sstemp.str(""); sstemp.clear(); sstemp<<setprecision(PRECISION)<<c.y;
|
||||
sstemp>>mid_id;
|
||||
vert_id = vert_id + " " + mid_id;
|
||||
sstemp.str(""); sstemp.clear(); sstemp<<setprecision(PRECISION)<<c.z;
|
||||
sstemp>>mid_id;
|
||||
vert_id = vert_id + " " + mid_id;
|
||||
return vert_id;
|
||||
}
|
||||
//测试打开输入文件 如果成功则返回0并输出信息 否则返回1
|
||||
int open_infile(ifstream &infile,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,char* filename){
|
||||
outfile.open(filename);
|
||||
if (!outfile){
|
||||
cerr << BOLDRED << "error ==> " << RESET << "fail to create the file: " << filename << endl;
|
||||
return -1;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
//规则网络插值 长方形内数据插值 距离平方反比
|
||||
/*长方体示意图*/
|
||||
// y
|
||||
// |
|
||||
// |
|
||||
// 3------------2
|
||||
// | |
|
||||
// | |
|
||||
// | |
|
||||
// 0------------1--->x
|
||||
// 左下角坐标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)
|
||||
{
|
||||
double res = 0;
|
||||
double total_dist = 0;
|
||||
double dist[4] = {0,0,0,0};
|
||||
double val[4];
|
||||
val[0] = d0; val[1] = d1; val[2] = d2; val[3] = d3;
|
||||
dist[0] = 1.0/(1e-30+(x-x0)*(x-x0)+(y-y0)*(y-y0));
|
||||
dist[1] = 1.0/(1e-30+(x-dx-x0)*(x-dx-x0)+(y-y0)*(y-y0));
|
||||
dist[2] = 1.0/(1e-30+(x-dx-x0)*(x-dx-x0)+(y-dy-y0)*(y-dy-y0));
|
||||
dist[3] = 1.0/(1e-30+(x-x0)*(x-x0)+(y-dy-y0)*(y-dy-y0));
|
||||
for (int i = 0; i < 4; i++){
|
||||
total_dist += dist[i];
|
||||
}
|
||||
for (int i = 0; i < 4; i++){
|
||||
res += val[i]*dist[i]/total_dist;
|
||||
}
|
||||
return res;
|
||||
}
|
||||
@ -1,93 +0,0 @@
|
||||
#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 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; ///< 整型到整形的映射
|
||||
//模型块体参数
|
||||
struct modelist{
|
||||
char mod_type[1024];
|
||||
char val_type[1024];
|
||||
char mod_para[1024];
|
||||
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); //规则网络插值
|
||||
#endif
|
||||
@ -1,6 +1,6 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::InitObs(char* obs_para){
|
||||
int GM3D::InitObs(const char* obs_para){
|
||||
obspoint temp_obs;
|
||||
string temp_str;
|
||||
stringstream temp_ss;
|
||||
@ -10,8 +10,8 @@ int GM3D::InitObs(char* obs_para){
|
||||
|
||||
//按格式解析参数 初始化观测位置 用于正演计算
|
||||
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);
|
||||
xmin = GCTL_MIN(xs,xe); xmax = GCTL_MAX(xs,xe);
|
||||
ymin = GCTL_MIN(ys,ye); ymax = GCTL_MAX(ys,ye);
|
||||
|
||||
x = xs;
|
||||
while(x >= xmin && x <= xmax){
|
||||
@ -29,7 +29,7 @@ int GM3D::InitObs(char* obs_para){
|
||||
//解析失败 按文件读入 用于反演使用或者正演计算
|
||||
else{
|
||||
ifstream infile;
|
||||
if (open_infile(infile,obs_para)) return -1;
|
||||
gctl::open_infile(infile,obs_para);
|
||||
|
||||
while(getline(infile,temp_str)){
|
||||
if (*(temp_str.begin()) == '#') continue;
|
||||
@ -41,19 +41,14 @@ int GM3D::InitObs(char* obs_para){
|
||||
temp_obs.val = temp_obs.dev = 0.0;
|
||||
obs_p_.push_back(temp_obs);
|
||||
}
|
||||
else{
|
||||
cerr << BOLDYELLOW << "ignored ==> " << RESET << "wrong input: " << temp_str << endl;
|
||||
continue;
|
||||
}
|
||||
else throw GCTL_WRONG_FILE_FORMAT;
|
||||
}
|
||||
}
|
||||
infile.close();
|
||||
}
|
||||
|
||||
if (obs_p_.empty()){
|
||||
cerr << BOLDRED << "error ==> " << RESET << "fail to initial observations with the parameter: " << obs_para << endl;
|
||||
return -1;
|
||||
}
|
||||
if (obs_p_.empty())
|
||||
throw GM3D_INVALID_OBS_PARAMETER;
|
||||
else{
|
||||
obs_num_ = obs_p_.size();
|
||||
}
|
||||
|
||||
296
src/main.cpp
296
src/main.cpp
@ -1,196 +1,118 @@
|
||||
#include "gm3d.h"
|
||||
#include "disp_help.h"
|
||||
|
||||
void display_help_info(char* program_name){
|
||||
string ex_name = program_name;
|
||||
string ex_usage = program_name;
|
||||
ex_usage += " {-o<output-msh-file> -m<model-list-file> \
|
||||
[-d<xs>/<dx>/<xe>/<ys>/<dy>/<ye>/<zs>/<dz>/<ze>] [-r]} \
|
||||
| {-i<input-msh-file> -p<observation-file>|<xs>/<dx>/<xe>/<ys>/<dy>/<ye>/<elevation> \
|
||||
-f<output-data-file> -tVz|Vzx|Vzy|Vzz|DeltaT|DeltaTx|DeltaTy|DeltaTz|Hax|Hay|Za \
|
||||
[-v<I0>/<D0>/<I>/<D>] [-n<noise-mean>/<noise-dev>]} -e<element-data-name> [-h]";
|
||||
int main(int argc, char* argv[])
|
||||
{
|
||||
gctl::flags_parser fp;
|
||||
fp.set_proname(argv[0]);
|
||||
fp.set_proinfo("3D model construction, forward and invert modeling of gravity and magnetic data under the Cartesian coordinates.");
|
||||
|
||||
DispHelp dh;
|
||||
dh.changeLayerOut(0,10);
|
||||
dh.addHeadInfo(ex_name,"1.0","3D forward modeling of gravity and magnetic data under the Cartesian coordinates.","Dr. Yi Zhang (zhangyiss@icloud.com). School of Earth Sciences, Zhejiang University");
|
||||
dh.addUsage(ex_usage);
|
||||
dh.addOption("Filename of the input Gmsh(.msh) model file for forward calculation.","-i");
|
||||
dh.addOption("Filename of the output Gmsh(.msh) model file built with given parameters.","-o");
|
||||
dh.addOption("Filename of the output observation file of gravity or magnetic data.","-f");
|
||||
dh.addOption("3D dimensions of the model space. the suffix 's' means the starting coordinate and 'e' represents the ending coordinate in axial directions. \
|
||||
'dx', 'dy' and 'dz' are step lengths. The default value is 10/20/990/10/20/990/10/20/490. \
|
||||
The axial orientation adopted by the program is a right-hand Cartesian system with the z-axis point vertical downward.","-d");
|
||||
dh.addOption("Model file that contains different types of model parameter. See instructions for formats of different model types. Note that any line starts with '#' will be skipped. \
|
||||
Four different types of models are supported by the program. Keyword parameters for different model types are shown as:","-m");
|
||||
dh.addOptionSec("1. regular_block add|replace|erase <physical_value> <xmin>/<xmax>/<ymin>/<ymax>/<zmin>/<zmax>");
|
||||
dh.addOptionSec("2. tilted_block add|replace|erase <physical_value> <xmin_z>/<xmax_z>/<ymin_z>/<ymax_z>/<zmin>/<xmin_Z>/<xmax_Z>/<ymin_Z>/<ymax_Z>/<zmax>");
|
||||
dh.addOptionSec("3. sphere add|replace|erase <physical_value> <xcen>/<ycen>/<zcen>/<xradius>/<yradius>/<zradius>");
|
||||
dh.addOptionSec("4. interface add|replace|erase/top|bot <physical_value> <filename>");
|
||||
dh.addOption("Element data name of the input/output Gmsh(.msh) file. Note that the name must be around by \"\".","-e");
|
||||
dh.addOption("Observation locations. You can either initialize the observation points from parameters or a file. \
|
||||
Each line of the file contain coordinates y(easting), x(northing) and z(elevation) of an observation point.","-p");
|
||||
dh.addOption("Forward component Vz, Vzx, Vzy or Vzz for gravitational data and \
|
||||
DeltaT, DeltaTx, DeltaTy, DeltaTz, Hax, Hay and Za for magnetic data.","-t");
|
||||
dh.addOption("Inclination and declination of the geomagnetic field and magnetization.","-v");
|
||||
dh.addOption("Add noise to the forward calculated data","-n");
|
||||
dh.addOption("Remove model elements with no data in the output Gmsh(.msh) file.","-r");
|
||||
dh.addOption("Display help information.","-h");
|
||||
dh.show();
|
||||
return;
|
||||
}
|
||||
fp.add_opt('i', "input-mesh", required_argument, NULL, "Input Gmsh(.msh) model file for forward calculation.", "<file>", false);
|
||||
fp.add_opt('o', "output-mesh", required_argument, NULL, "Output Gmsh(.msh) model file for model construction.", "<file>", false);
|
||||
fp.add_opt('a', "anomalies", required_argument, NULL, "Output anomalies of the gravity or magnetic data.", "<file>", false);
|
||||
fp.add_opt('d', "dimension", required_argument, NULL, "3D dimensions of the model space. the suffix 's' means the starting coordinate \
|
||||
and 'e' represents the ending coordinate in axial directions. 'dx', 'dy' and 'dz' are step lengths. The default value is 10/20/990/10/20/990/10/20/490. \
|
||||
The axial orientation adopted by the program is a right-hand Cartesian system with the z-axis point vertical downward.",
|
||||
"<xs>/<dx>/<xe>/<ys>/<dy>/<ye>/<zs>/<dz>/<ze>", false);
|
||||
fp.add_opt('m', "model", required_argument, NULL, "Model file that contains model parameters. See instructions for formats of different model types. \
|
||||
Note that any line starts with '#' will be skipped. Four different types of models are supported by the program. Keyword parameters for different model types are shown as: \
|
||||
1. regular_block add|replace|erase <physical_value> <xmin>/<xmax>/<ymin>/<ymax>/<zmin>/<zmax> 2. tilted_block add|replace|erase <physical_value> <xmin_z>/<xmax_z>/<ymin_z>/<ymax_z>/<zmin>/<xmin_Z>/<xmax_Z>/<ymin_Z>/<ymax_Z>/<zmax> \
|
||||
3. sphere add|replace|erase <physical_value> <xcen>/<ycen>/<zcen>/<xradius>/<yradius>/<zradius> 4. interface add|replace|erase/top|bot <physical_value> <filename>",
|
||||
"<file>", false);
|
||||
fp.add_opt('e', "element", required_argument, NULL, "Element data name of the input/output Gmsh(.msh) file. Note that the name must be around by \"\".",
|
||||
"<name>", false);
|
||||
fp.add_opt('p', "position", required_argument, NULL, "Observation positions. You can either initialize the observation points from parameters or a file. \
|
||||
Each line of the file contain coordinates y(easting), x(northing) and z(elevation) of an observation point.",
|
||||
"<observation-file>|<xs>/<dx>/<xe>/<ys>/<dy>/<ye>/<elevation>", false);
|
||||
fp.add_opt('t', "type", required_argument, NULL, "Forward component Vz, Vzx, Vzy or Vzz for gravitational data and \
|
||||
DeltaT, DeltaTx, DeltaTy, DeltaTz, Hax, Hay and Za for magnetic data.",
|
||||
"Vz|Vzx|Vzy|Vzz|DT|DTx|DTy|DTz|Hax|Hay|Za", false);
|
||||
fp.add_opt('v', "vector", required_argument, NULL, "Inclination and declination of the geomagnetic field and magnetization.", "<I0>/<D0>/<I>/<D>", false);
|
||||
fp.add_opt('n', "noise", required_argument, NULL, "Add noise to the forward calculated data", "<mean>/<std>", false);
|
||||
fp.add_opt('r', "remove", no_argument, NULL, "Remove model elements with no data in the output Gmsh(.msh) file.", 0, false);
|
||||
fp.add_opt('l', "log", required_argument, NULL, "Log file. A default name will be used if this option is not set.", "<file>", false);
|
||||
fp.add_opt('h', "help", no_argument, NULL, "Show help information.", 0, false);
|
||||
fp.configure(argc, argv);
|
||||
|
||||
int main(int argc, char* argv[]){
|
||||
GM3D gm3d_instance;
|
||||
char dimension[1024] = "10/20/990/10/20/990/10/20/490";
|
||||
char modelname[1024] = "NULL";
|
||||
char in_mshname[1024] = "NULL";
|
||||
char out_mshname[1024] = "Untitled.msh";
|
||||
char elename[1024] = "Untitled";
|
||||
char obsername[1024] = "NULL";
|
||||
char out_obserfile[1024] = "NULL";
|
||||
char noise_para[1024] = "NULL";
|
||||
char forward_type[1024] = "Vz";
|
||||
char mag_field_para[1024] = "90/0/90/0";
|
||||
bool build_model = true;
|
||||
bool remove_null = false;
|
||||
|
||||
opterr = 0; //内置参数 若不为0则会在发生遭遇错误时输出一条信息到屏幕
|
||||
|
||||
if (argc == 1){
|
||||
display_help_info(argv[0]);
|
||||
return 0;
|
||||
}
|
||||
|
||||
int curr;
|
||||
/*循环拾取参数 最后一个参数为-1 需要变量的参数后跟一个冒号 可有可无参数跟两个冒号*/
|
||||
while((curr = getopt(argc,argv,"hrd:m:i:o:e:p:f:n:t:v:")) != -1){
|
||||
/*匹配命令*/
|
||||
switch (curr){
|
||||
case 'h': //显示帮助信息
|
||||
display_help_info(argv[0]);
|
||||
return 0;
|
||||
case 'r':
|
||||
remove_null = true;
|
||||
break;
|
||||
case 'd':
|
||||
if (1!=sscanf(optarg,"%s",dimension)){
|
||||
cout << "error ==> wrong format of " << optarg << endl;
|
||||
}
|
||||
break;
|
||||
case 'm':
|
||||
if (1!=sscanf(optarg,"%s",modelname)){
|
||||
cout << "error ==> wrong format of " << optarg << endl;
|
||||
}
|
||||
break;
|
||||
case 'i':
|
||||
if (1!=sscanf(optarg,"%s",in_mshname)){
|
||||
cout << "error ==> wrong format of " << optarg << endl;
|
||||
}
|
||||
else build_model = false;
|
||||
break;
|
||||
case 'o':
|
||||
if (1!=sscanf(optarg,"%s",out_mshname)){
|
||||
cout << "error ==> wrong format of " << optarg << endl;
|
||||
}
|
||||
break;
|
||||
case 'f':
|
||||
if (1!=sscanf(optarg,"%s",out_obserfile)){
|
||||
cout << "error ==> wrong format of " << optarg << endl;
|
||||
}
|
||||
break;
|
||||
case 'n':
|
||||
if (1!=sscanf(optarg,"%s",noise_para)){
|
||||
cout << "error ==> wrong format of " << optarg << endl;
|
||||
}
|
||||
break;
|
||||
case 'e':
|
||||
if (1!=sscanf(optarg,"%s",elename)){
|
||||
cout << "error ==> wrong format of " << optarg << endl;
|
||||
}
|
||||
break;
|
||||
case 'v':
|
||||
if (1!=sscanf(optarg,"%s",mag_field_para)){
|
||||
cout << "error ==> wrong format of " << optarg << endl;
|
||||
}
|
||||
break;
|
||||
case 'p':
|
||||
if (1!=sscanf(optarg,"%s",obsername)){
|
||||
cout << "error ==> wrong format of " << optarg << endl;
|
||||
}
|
||||
break;
|
||||
case 't':
|
||||
if (1!=sscanf(optarg,"%s",forward_type)){
|
||||
cout << "error ==> wrong format of " << optarg << endl;
|
||||
}
|
||||
break;
|
||||
case '?': //处理未定义或错误参数
|
||||
if (optopt == 'd' || optopt == 'm' || optopt == 'o'
|
||||
|| optopt == 'e' || optopt == 'i' || optopt == 'p'
|
||||
|| optopt == 'n' || optopt == 't' || optopt == 'v'){
|
||||
fprintf (stderr, "Option -%c requires an argument.\n", optopt);
|
||||
return -1;
|
||||
}
|
||||
else if (isprint(optopt)){
|
||||
fprintf (stderr, "Unknown option `-%c'.\n", optopt);
|
||||
return -1;
|
||||
}
|
||||
else{
|
||||
fprintf (stderr,"Unknown option character `\\x%x'.\n",optopt);
|
||||
return -1;
|
||||
}
|
||||
break;
|
||||
default:
|
||||
abort();
|
||||
}
|
||||
}
|
||||
|
||||
if (build_model){
|
||||
if (gm3d_instance.BuildRegularGrid(dimension)) return 0;
|
||||
if (gm3d_instance.AddModels(modelname)) return 0;
|
||||
gm3d_instance.RegisteredOuput(remove_null);
|
||||
if (gm3d_instance.OutMshFile(out_mshname,elename)) return 0;
|
||||
}
|
||||
else{
|
||||
if (gm3d_instance.ReadModel(in_mshname,elename)) return 0;
|
||||
if (gm3d_instance.InitObs(obsername)) return 0;
|
||||
if (!strcmp(forward_type,"Vz")){
|
||||
gm3d_instance.ForwardVz(noise_para);
|
||||
}
|
||||
else if (!strcmp(forward_type,"Vzx")){
|
||||
gm3d_instance.ForwardVzx(noise_para);
|
||||
}
|
||||
else if (!strcmp(forward_type,"Vzy")){
|
||||
gm3d_instance.ForwardVzy(noise_para);
|
||||
}
|
||||
else if (!strcmp(forward_type,"Vzz")){
|
||||
gm3d_instance.ForwardVzz(noise_para);
|
||||
}
|
||||
else if (!strcmp(forward_type,"DeltaT")){
|
||||
gm3d_instance.ForwardDeltaT(noise_para,mag_field_para);
|
||||
}
|
||||
else if (!strcmp(forward_type,"DeltaTx")){
|
||||
gm3d_instance.ForwardDeltaTx(noise_para,mag_field_para);
|
||||
}
|
||||
else if (!strcmp(forward_type,"DeltaTy")){
|
||||
gm3d_instance.ForwardDeltaTy(noise_para,mag_field_para);
|
||||
}
|
||||
else if (!strcmp(forward_type,"DeltaTz")){
|
||||
gm3d_instance.ForwardDeltaTz(noise_para,mag_field_para);
|
||||
}
|
||||
else if (!strcmp(forward_type,"Hax")){
|
||||
gm3d_instance.ForwardHax(noise_para,mag_field_para);
|
||||
}
|
||||
else if (!strcmp(forward_type,"Hay")){
|
||||
gm3d_instance.ForwardHay(noise_para,mag_field_para);
|
||||
}
|
||||
else if (!strcmp(forward_type,"Za")){
|
||||
gm3d_instance.ForwardZa(noise_para,mag_field_para);
|
||||
}
|
||||
else{
|
||||
cerr << BOLDYELLOW << "error ==> " << RESET << "wrong forward component type: " << forward_type << endl;
|
||||
try
|
||||
{
|
||||
if(argc == 1 || fp.set_opt('h'))
|
||||
{
|
||||
fp.show_help_page();
|
||||
return 0;
|
||||
}
|
||||
if (gm3d_instance.OutObs(out_obserfile)) return 0;
|
||||
|
||||
bool build_model = true;
|
||||
bool remove_null = false;
|
||||
std::string modelname, in_mshname, out_mshname, elename, obsername, out_obserfile;
|
||||
std::string noise_para, forward_type, mag_field_para, dimension, logfile_name;
|
||||
|
||||
fp.get_argv({'m', 'i', 'o', 'e', 'p', 'a', 'n', 't', 'v', 'd', 'l'},
|
||||
{&modelname, &in_mshname, &out_mshname, &elename, &obsername, &out_obserfile,
|
||||
&noise_para, &forward_type, &mag_field_para, &dimension, &logfile_name});
|
||||
if (logfile_name == "NULL") logfile_name = "gm3d.history.txt";
|
||||
if (in_mshname != "NULL") build_model = false;
|
||||
if (!fp.set_opt('r')) remove_null = true;
|
||||
|
||||
// 查看是否通过强制参数检查
|
||||
if (!fp.pass_mandatory()) return 0;
|
||||
|
||||
GM3D gm3d_instance;
|
||||
if (build_model)
|
||||
{
|
||||
gm3d_instance.BuildRegularGrid(dimension.c_str());
|
||||
gm3d_instance.AddModels(modelname.c_str());
|
||||
gm3d_instance.RegisteredOuput(remove_null);
|
||||
gm3d_instance.OutMshFile(out_mshname.c_str(), elename);
|
||||
}
|
||||
else
|
||||
{
|
||||
gm3d_instance.ReadModel(in_mshname.c_str(),elename.c_str());
|
||||
gm3d_instance.InitObs(obsername.c_str());
|
||||
|
||||
if (forward_type == "Vz") gm3d_instance.ForwardVz(noise_para.c_str());
|
||||
else if (forward_type == "Vzx") gm3d_instance.ForwardVzx(noise_para.c_str());
|
||||
else if (forward_type == "Vzy") gm3d_instance.ForwardVzy(noise_para.c_str());
|
||||
else if (forward_type == "Vzz") gm3d_instance.ForwardVzz(noise_para.c_str());
|
||||
else if (forward_type == "DT") gm3d_instance.ForwardDeltaT(noise_para.c_str(),mag_field_para.c_str());
|
||||
else if (forward_type == "DTx") gm3d_instance.ForwardDeltaTx(noise_para.c_str(),mag_field_para.c_str());
|
||||
else if (forward_type == "DTy") gm3d_instance.ForwardDeltaTy(noise_para.c_str(),mag_field_para.c_str());
|
||||
else if (forward_type == "DTz") gm3d_instance.ForwardDeltaTz(noise_para.c_str(),mag_field_para.c_str());
|
||||
else if (forward_type == "Hax") gm3d_instance.ForwardHax(noise_para.c_str(),mag_field_para.c_str());
|
||||
else if (forward_type == "Hay") gm3d_instance.ForwardHay(noise_para.c_str(),mag_field_para.c_str());
|
||||
else if (forward_type == "Za") gm3d_instance.ForwardZa(noise_para.c_str(),mag_field_para.c_str());
|
||||
else throw "wrong forward component type.";
|
||||
|
||||
gm3d_instance.OutObs(out_obserfile.c_str());
|
||||
}
|
||||
|
||||
//最后输出一个记录文件 我们使用唯一的输出文件名称 所以每次在统一目录运行都会覆盖
|
||||
ofstream log_out;
|
||||
gctl::open_outfile(log_out, logfile_name);
|
||||
|
||||
time_t now = time(0);
|
||||
char* dt = ctime(&now);
|
||||
|
||||
log_out << "Time: " << dt;
|
||||
log_out << "Mesh: " << dimension << endl;
|
||||
log_out << "Model id: " << gm3d_instance.model_id_ << endl;
|
||||
log_out << "Input .msh file: " << in_mshname << endl;
|
||||
log_out << "Output .msh file: " << out_mshname << endl;
|
||||
log_out << "Model list file: " << modelname << endl;
|
||||
log_out << "Model element name: " << elename << endl;
|
||||
log_out << "observation parameters: " << obsername << endl;
|
||||
log_out << "Output observation file: " << out_obserfile << endl;
|
||||
log_out << "noise level: " << noise_para << endl;
|
||||
log_out << "forward type: " << forward_type << endl;
|
||||
log_out << "magnetic field parameters: " << mag_field_para << endl;
|
||||
log_out.close();
|
||||
}
|
||||
catch(std::string err_str)
|
||||
{
|
||||
GCTL_ShowWhatError(err_str, GCTL_ERROR_ERROR, 0, 0, 0);
|
||||
}
|
||||
catch(const char* err_char)
|
||||
{
|
||||
GCTL_ShowWhatError(err_char, GCTL_ERROR_ERROR, 0, 0, 0);
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
@ -1,10 +1,22 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::OutMshFile(char* filename,string data_name){
|
||||
int GM3D::OutMshFile(const char* filename,string data_name){
|
||||
time_t now = time(0);
|
||||
char* dt = ctime(&now);
|
||||
|
||||
if (!strcmp(filename,"NULL")) return 0;
|
||||
|
||||
ofstream outfile;
|
||||
if (open_outfile(outfile,filename)) return -1;
|
||||
gctl::open_outfile(outfile,filename);
|
||||
|
||||
gctl::random_char(12, model_id_);
|
||||
// 输出一个注释信息单元 包含建模的时间 与 一个唯一的十二位编码
|
||||
// 这个编码会同时出现在记录文件与正演结果文件中 通过辨别这些文件中的编码可确定其对应关系
|
||||
outfile << "$Comments" << endl;
|
||||
outfile << "Created by : gm3d" << endl;
|
||||
outfile << "At time: " << dt; //dt包含换行符
|
||||
outfile << "Model id: " << model_id_ << endl;
|
||||
outfile << "$EndComments" << endl;
|
||||
|
||||
//好啦 我们这里输出的模型类型应该是块体
|
||||
outfile<<"$MeshFormat"<<endl<<"2.2 0 8"<<endl<<"$EndMeshFormat"<<endl<<"$Nodes"<<endl<< out_vert_ids_.size() <<endl;
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::OutNeighborFile(char* v_name,char* b_name){
|
||||
int GM3D::OutNeighborFile(const char* v_name,const char* b_name){
|
||||
/*
|
||||
if (strcmp(v_name,"NULL") || strcmp(b_name,"NULL")){
|
||||
//整理块体间的相邻关系 先初始化顶点相邻数组
|
||||
|
||||
@ -1,8 +1,15 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::OutObs(char* filename){
|
||||
int GM3D::OutObs(const char* filename){
|
||||
time_t now = time(0);
|
||||
char* dt = ctime(&now);
|
||||
|
||||
ofstream outfile;
|
||||
if (open_outfile(outfile,filename)) return -1;
|
||||
gctl::open_outfile(outfile,filename);
|
||||
|
||||
outfile << "# Created by : gm3d" << endl;
|
||||
outfile << "# At time: " << dt; //dt包含换行符
|
||||
outfile << "# Model id: " << model_id_ << endl;
|
||||
outfile << "# y(m)-easting x(m)-northing ele(m) obs-val(mGal|Eo) obs-dev(mGal|Eo)" << endl;
|
||||
for (int i = 0; i < obs_num_; i++){
|
||||
outfile << obs_p_[i].y << " " << obs_p_[i].x << " " << -1.0*obs_p_[i].z << " "
|
||||
|
||||
@ -1,114 +0,0 @@
|
||||
//#ifdef _WINDOWS
|
||||
//#include <windows.h>
|
||||
//#else
|
||||
//#include <sys/ioctl.h>
|
||||
//#endif
|
||||
|
||||
#include "progress_bar.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;
|
||||
}
|
||||
|
||||
}
|
||||
@ -1,41 +0,0 @@
|
||||
#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
|
||||
#define CHARACTER_WIDTH_PERCENTAGE 4
|
||||
|
||||
class ProgressBar
|
||||
{
|
||||
public:
|
||||
|
||||
ProgressBar();
|
||||
ProgressBar(unsigned long n_, const char *description_="", std::ostream& out_=std::cerr);
|
||||
|
||||
void SetFrequencyUpdate(unsigned long frequency_update_);
|
||||
void SetStyle(const char* unit_bar_, const char* unit_space_);
|
||||
|
||||
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
|
||||
@ -1,26 +1,27 @@
|
||||
#include "gm3d.h"
|
||||
|
||||
int GM3D::ReadModel(char* filename,char* input_forward_model_name){
|
||||
int GM3D::ReadModel(const char* filename,const char* input_forward_model_name){
|
||||
int temp_int,ele_type,attri_num,temp_attri,temp_id;
|
||||
double temp_val;
|
||||
_1dArray temp_model;
|
||||
cpoint temp_vert;
|
||||
gctl::_1d_vector temp_model;
|
||||
vertex temp_vert;
|
||||
cube temp_cu;
|
||||
string temp_str;
|
||||
stringstream temp_ss;
|
||||
|
||||
ifstream mshin;
|
||||
if (open_infile(mshin,filename)) return -1; //检查并打开模型文件
|
||||
gctl::open_infile(mshin,filename);
|
||||
|
||||
while(getline(mshin,temp_str)){
|
||||
//读入模型空间顶点集 msh文件版本为2.2
|
||||
if (temp_str == "$Nodes"){
|
||||
getline(mshin,temp_str);
|
||||
temp_ss = str2ss(temp_str);
|
||||
gctl::str2ss(temp_str, temp_ss);
|
||||
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);
|
||||
gctl::str2ss(temp_str, temp_ss);
|
||||
temp_ss >> temp_vert.id >> temp_vert.x >> temp_vert.y >> temp_vert.z;
|
||||
model_vert_[i] = temp_vert;
|
||||
}
|
||||
@ -28,12 +29,12 @@ int GM3D::ReadModel(char* filename,char* input_forward_model_name){
|
||||
//读入模型空间单元体
|
||||
else if (temp_str == "$Elements"){
|
||||
getline(mshin,temp_str);
|
||||
temp_ss = str2ss(temp_str);
|
||||
gctl::str2ss(temp_str, temp_ss);
|
||||
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);
|
||||
gctl::str2ss(temp_str, temp_ss);
|
||||
temp_ss >> temp_cu.cen.id >> ele_type;
|
||||
//只读入块体
|
||||
if (ele_type == 5){
|
||||
@ -46,12 +47,22 @@ int GM3D::ReadModel(char* filename,char* input_forward_model_name){
|
||||
}
|
||||
}
|
||||
}
|
||||
//读入model_id
|
||||
else if (temp_str == "$Comments"){
|
||||
while (strcmp(temp_str.c_str(),"$EndComments")){
|
||||
getline(mshin,temp_str);
|
||||
if ( 1 == sscanf(temp_str.c_str(),"Model id: %s",model_id_)){
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
else continue; //不能识别的单元都被忽略了
|
||||
}
|
||||
mshin.close();
|
||||
|
||||
//第二次读入模型文件 初始化模型单元属性
|
||||
if (open_infile(mshin,filename)) return -1; //检查并打开模型文件
|
||||
gctl::open_infile(mshin,filename);
|
||||
|
||||
while(getline(mshin,temp_str)){
|
||||
//读入模型单元属性 注意因为msh文件中$ElementData并未注明所属元素类型
|
||||
//所以可能会将其他元素类型的属性值也读入 但因为其在pyIdMap中并未注册 所以属性值会全为0 在后续使用时我们需要通过名称辨别
|
||||
@ -62,11 +73,11 @@ int GM3D::ReadModel(char* filename,char* input_forward_model_name){
|
||||
input_model_names_.push_back(temp_str);
|
||||
for (int i = 0; i < 6; i++) //跳过元素属性前面的值 最后一次为当前元素块的个数
|
||||
getline(mshin,temp_str);
|
||||
temp_ss = str2ss(temp_str);
|
||||
gctl::str2ss(temp_str, temp_ss);
|
||||
temp_ss >> temp_int;
|
||||
for (int i = 0; i < temp_int; i++){
|
||||
getline(mshin,temp_str);
|
||||
temp_ss = str2ss(temp_str);
|
||||
gctl::str2ss(temp_str, temp_ss);
|
||||
temp_ss >> temp_id >> temp_val; //读入单元体索引与属性值
|
||||
temp_model[temp_id] = temp_val;
|
||||
}
|
||||
@ -87,15 +98,23 @@ int GM3D::ReadModel(char* filename,char* input_forward_model_name){
|
||||
char forward_model_name[1024] = "\"";
|
||||
strcat(forward_model_name,input_forward_model_name);
|
||||
strcat(forward_model_name,"\"");
|
||||
for (int i = 0; i < input_model_names_.size(); i++){
|
||||
if (!strcmp(input_model_names_[i].c_str(),forward_model_name)){
|
||||
forward_model_ = input_models_[i];
|
||||
clog << "forward model initialized from file." << endl;
|
||||
|
||||
try
|
||||
{
|
||||
for (int i = 0; i < input_model_names_.size(); i++){
|
||||
if (!strcmp(input_model_names_[i].c_str(),forward_model_name)){
|
||||
forward_model_ = input_models_[i];
|
||||
throw GM3D_FORWARD_MODEL_INITIALIZED;
|
||||
}
|
||||
}
|
||||
}
|
||||
catch(const char* err_str)
|
||||
{
|
||||
GCTL_ShowWhatError(err_str, GCTL_MESSAGE_ERROR, 0, 0, 0);
|
||||
}
|
||||
|
||||
//计算块体的中心位置和尺寸
|
||||
cpoint corner[8];
|
||||
vertex corner[8];
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
for (int j = 0; j < 8; j++){
|
||||
corner[j] = model_vert_[model_cube_[i].ids[j]];
|
||||
|
||||
@ -7,7 +7,7 @@ int GM3D::RegisteredOuput(bool remove_empty_element){
|
||||
count = 0;
|
||||
//遍历所有块体数据 注册有值的块体
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
if (model_block_val_[i] != BDL_MAX){
|
||||
if (model_block_val_[i] != GCTL_BDL_MAX){
|
||||
out_ele_data_ids_.push_back(i);
|
||||
out_ele_ids_.push_back(i);
|
||||
ele_data_out_map_[count] = count;
|
||||
@ -49,7 +49,7 @@ int GM3D::RegisteredOuput(bool remove_empty_element){
|
||||
//注册所有有值的块体数据
|
||||
count = 0;
|
||||
for (int i = 0; i < model_num_; i++){
|
||||
if (model_block_val_[i] != BDL_MAX){
|
||||
if (model_block_val_[i] != GCTL_BDL_MAX){
|
||||
ele_data_out_map_[count] = i;
|
||||
out_ele_data_ids_.push_back(i);
|
||||
count++;
|
||||
|
||||
Loading…
Reference in New Issue
Block a user