大地电磁正演程序MT2D主程序分析

0 相关背景信息

见大地电磁二维正演程序--详细介绍（https://blog.csdn.net/spvfly/article/details/92795423）。

程序的框架如下图

1 主程序MT2D.cpp代码分析

01 头文件部分

包含四个建好的类

#include "mesh2d.h"   //网格处理类  包含读入网格数据等功能
#include "dofs.h"         //有限单元的自由度计算
#include "2dmt.h"        //有限单元的计算    如 Ke = p的大型稀疏矩阵  计算出全局的电场（e）和磁场(h)
#include "post.h"         //后处理    提取地面测量点出的  电场（e）或  磁场（h）的值以及相位  以便利用matlab可视化作图

02 全局变量/函数

/*

全局变量region_table存储 人为划分的 网格区域属性（config文件中 设置，mesh2d.h类中有相应的读入方法），其格式如下：

config文件中的区域信息格式如下
2D-1.1  网格文件名（无后缀的）
6      有6个区域

标号    σ     ε    μ
 -1   1E-10  1    1
 1    0.04     1    1
 2     0.1      1    1
 3     0.4      1    1
 4    0.001   1    1
 5     0.2      1    1  

*/

std::map<int, std::vector<double> > region_table;  

// 函数electric_parameters，用于调取 三角单元 的 物理属性（电性）参数

void electric_parameters(double& cond, double& epsilon, double& mu,
                         int marker)
{
  assert(!region_table.empty());
  typedef std::map<int, std::vector<double> >::iterator IT;
  IT it= region_table.find(marker);
  assert(it!=region_table.end());
  
  std::vector<double>& parameter = (*it).second;
  assert(parameter.size()==3);  
  
  cond    = parameter[0];    
  epsilon = parameter[1]*EM::epsilon0;
  mu      = parameter[2]*EM::mu0;
  
  return;
}

 

//主程序入口

int main(int argc, char *argv[]) {
  
  /* 
     0, start model   ==》读入模型信息 在xx.config文件中
  */  
  if (argc <2 ) {
    printf("Usage: %s input_paramters_filename\n",argv[0]);
    return 1;    
  }  
  const char *filename = argv[1];  
  std::ifstream in_stream (filename);
  assert(in_stream.good());
  
  std::string starting_model;
  double frequency, theta_in;
  unsigned int algorithm;
  int    station_marker;
  unsigned int n_layers;
  std::vector<std::vector<double> > EP;
  unsigned int n_regions;

  // Part 0
  in_stream >> starting_model;
  std::cout<<starting_model<<"\n";
  // Part I
  in_stream >> frequency 
            >> theta_in
            >> algorithm
            >> station_marker;
  std::cout<<frequency<<"\n"
           <<theta_in<<"\n"
           << algorithm<<"\n"
           << station_marker<<"\n";
  // Part II
  in_stream >> n_layers;
  std::cout<< n_layers<<"\n";
  EP.resize(n_layers);
  for(int i=0; i<n_layers; i++) {
    EP[i].resize(4);
    for(int j=0; j<4; j++) { 
      double temp;
      in_stream >> temp;
      std::cout<<temp<<"\t";
      EP[i][j] = temp;
    }
    std::cout<<"\n";
  }
  // Part III  
  in_stream >> n_regions;
  std::cout<<n_regions<<"\n";
  for(int i=0; i<n_regions; i++)  {
    std::vector<double> temp;
    int marker;
    in_stream >> marker;
    std::cout<<marker<<"\t";
    temp.resize(3);
    for(int j=0; j<3; j++) {
      in_stream >> temp[j];
      std::cout<<temp[j]<<"\t";
    }
    std::cout<<"\n";
    region_table[marker] = temp;
  }
  
  /* 
     1, generate mesh  ==》生成网格  读取triangle生成的网格信息
  */
  Mesh2D mesh(starting_model);
  mesh.libmesh_read_triangle_files();

  /* 
     2, generate boundary condition   ==》边界条件 
  */
  BC bc_te(EP, frequency, theta_in, EM::TE); // TE source
  BC bc_tm(EP, frequency, theta_in, EM::TM); // TM source

  /* 
     3, generate dofs class  ==》分析 每个单元的  自由度
  */
  DOFs dofs(mesh);
  dofs.distribute_dofs();   // using triangle's order

  /* 
     4, generate MT class to compute Ex and Hx  计算 电场 和 磁场
  */
  MT mt(argc, argv, mesh, dofs, frequency, bc_te, bc_tm, algorithm); // 
  mt.solve();
  DenseVector Ex, Hx;
  mt.get_X(Ex, Hx); // solution for TE and TM mode
  

  /* 
     5, output to matlab for visualization   后处理  用matlab成图
  */
  POST post(&mesh, &dofs, Ex, Hx, frequency);
  post.identify_mt_stations(station_marker);
  
  return 0;

}

注：本程序为任政勇教授开源的MT2D程序（https://sourceforge.net/projects/mt2d/），只是在学习其编程的思路，不涉及其他目的。