有限要素法(FEM)のデータモデル EXODUS導入の覚書

はじめに

EXODUSは有限要素解析のデータ格納用のデータファイルです。
他のコードで計算結果を流用する際の中間生成のファイルとしても、計算結果を格納する形式としても利用できるフォーマットになります。

これまで自作のFEMについては計算結果をvtk形式で出力していましたが、EXODUSの方が再利用性が高そうだったのでこちらに移行中です。

導入

下記リポジトリよりコードを持ってきます。
https://github.com/gsjaardema/seacas#exodus
導入についてもこちらに則れば基本的に問題ないかと思います。

ここでは~/localでチェックアウト・ビルドすることを想定します。

cd ~
mkdir local
cd local
git clone https://github.com/gsjaardema/seacas.git
cd seacas && export ACCESS=`pwd`

Third-Party Library (TPL)としてNetCDFとHDF5を用意する必要がありますので、それぞれ下記からダウンロードします。
https://www.unidata.ucar.edu/downloads/netcdf/
https://portal.hdfgroup.org/display/support/HDF5%201.12.0

cd TPL/hdf5
wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.12/hdf5-1.12.0/src/hdf5-1.12.0.tar.gz
tar -zxvf ./hdf5-1.12.0.tar.gz
cd hdf5-1.12.0
../runconfigure.sh
cd ../../netcdf
wget https://www.unidata.ucar.edu/downloads/netcdf/ftp/netcdf-c-4.7.4.tar.gz
tar -zxvf ./netcdf-c-4.7.4.tar.gz
cd netcdf-c-4.7.4/
mkdir build
cd build
sh ../../runcmake.sh
make && make install
cd ../../../../
mkdir build
../cmake-exodus
make && make install

使い方

基本的には下記マニュアルを見ればOK
https://gsjaardema.github.io/seacas/exodusII-new.pdf

サンプルとして、寸法が1 x 1 x 10の片持ち梁を10個の8節点6面体要素で分割したメッシュをEXODUS形式で出力する場合のコードを書きに示します。

#include <exodusII.h>
#include <iostream>
#include <vector>

int main(){
  int N;
  N = 10; // the length of cantilever

  int CPU_word_size, IO_word_size, exoid;
  CPU_word_size = sizeof(double);
  IO_word_size = 8;

  exoid = ex_create("test.exo",                //file Name
            EX_CLOBBER,                //create mode
            &CPU_word_size,            //CPU float word size in bytes
            &IO_word_size);            //I/O float word size in bytes


  int dim = 3;           // dimension
  int NoN = 4 + 4 * N;   // num of nodes
  int NoE = N;       // num of elements
  int NoEB = 1;          // num of element blocks
  int NoNS = 0;          // num of node sets
  int NoSS = 0;          // num of side sets

  ex_put_init(exoid, "Test", dim, NoN, NoE, NoEB, NoNS, NoSS);
  std::vector<double> x(NoN), y(NoN), z(NoN), i2nn(NoN);
  char* coord_names[3];
  coord_names[0] = "xcoor";
  coord_names[1] = "ycoor";
  coord_names[2] = "zcoor";

  for(int i=0;i<=N;i++){
    x[i*4] = 0;   y[i*4] = 0;   z[i*4]=i;   i2nn[i*4] = i*4 + 1;
    x[i*4+1] = 1; y[i*4+1] = 0; z[i*4+1]=i; i2nn[i*4+1] = i*4 + 2;
    x[i*4+2] = 1; y[i*4+2] = 1; z[i*4+2]=i; i2nn[i*4+2] = i*4 + 3;
    x[i*4+3] = 0; y[i*4+3] = 1; z[i*4+3]=i; i2nn[i*4+3] = i*4 + 4;
  }
  ex_put_coord(exoid, x.data(), y.data(), z.data());
  ex_put_coord_names(exoid, coord_names);
  ex_put_node_num_map(exoid, i2nn.data());


  //for(int i=0;i<NoE;i++)i2nn[i] = i+1;

  int NoN_for_EB = NoE; // number of nodes for element block
  int NoN_per_Elem = 8; // number of nodes of each element
  int ebid = 1;         // element block id
  ex_put_elem_block(exoid, ebid, "hex", NoN_for_EB, NoN_per_Elem, 0);

  std::vector<int> connectivity(NoE * NoN_per_Elem);
  std::vector<int> i2en(NoE); // index to element number

  for(int i=0; i<NoE; i++){
    for(int j=0;j<8;j++){
      connectivity[i*8 + j] = i*4 + j +1;
    }
    i2en[i] = i+1;
  }
  ex_put_elem_conn(exoid, ebid, connectivity.data());
  ex_put_elem_num_map(exoid, i2en.data());
  ex_close(exoid);

  return 0;
}

例えば下記のようにコンパイルできます。

g++ main.cpp -I ~/local/seacas/include -L ~/local/seacas/lib -lexodus

出力したtest.exoをParaviewで見ると下記のようになります。

今回は節点量や要素量を出力していませんので、コンターは要素IDになります。

終わりに

今回は有限要素法の結果を格納可能なEXODUS形式について、その導入と動作確認を示しました。
もっと優れたファイルタイプがありましたら教えていただけると幸いです。

以上です。