以下が、今回のプログラムで使う木構造のためのクラス定義である。
pt minus 1 pt
#ifndef BHTREE_H
# define BHTREE_H
/*-----------------------------------------------------------------------------
* BHtree : basic class for C++ implementation of BH treecode
* J. Makino 2001/1/23
*-----------------------------------------------------------------------------
*/
#include "vector.h"
#include "particle.h"
class bhnode
{
private:
vector cpos;
real l;
bhnode * child[8];
particle * pfirst;
int nparticle;
vector pos;
real mass;
public:
bhnode(){
cpos = 0.0;
l = 0.0;
for(int i = 0; i<8;i++)child[i] = NULL;
pfirst = NULL;
nparticle = 0;
pos = 0.0;
mass = 0.0;
}
void create_tree_recursive(bhnode * & heap_top, int & heap_remainder);
void assign_root(vector root_pos, real length, particle * p, int nparticle);
void dump(int indent = 0);
int sanity_check();
void set_cm_quantities();
void accumulate_force_from_tree(vector & ipos, real eps2, real theta2,
vector & acc,
real & phi);
};
#endif
BH というのがやたらでてくるが、これはこのアルゴリズムを1986 年に提唱したBarnes
と Hut （当時はどちらもプリンストン高等研究所。今は Barnes はハワイ大
学天文学科）の頭文字をとっているからである。 cpos は立方体の中心
座標、 l は辺の長さである。 child は子セル（ノード）へのポ
インタの配列であり、最大 8 個の子セルがあるので要素数 8 の配列にしてあ
る。これを使うことで、線形リストのように1つだけを指すのではなく、8個の
場所を指せるわけで、それによって木構造が表現できるわけである。
pfirst は中にある粒子のリストの先頭を指す。 nparticle はこ のリストの長さである。
ここで、ちょっと注意して欲しいのは、長さは本来必要ではないということで ある。つまり、リストを作る時に、最後の要素が NULL を指していると いうふうに決めたので、長さはリストをたどっていって数えればわかるわけで、 本来長さを別に覚えておく必要はない。が、数えなくても長さがわかると便利 なことが多いので別に覚えておくことにする。このようにした場合には、情報 が冗長になっているので2つの情報が矛盾しないように、あるいはちゃんと意 味があるほうを使うように注意する必要がある。
実際に木を作る手順を見てみる。最初にすることは、1つ立方体をつくって、 それをルート (根）とし、全粒子をその中にあるということにする、具体的に は、全粒子をルートのpfirst から始まるリストに入れておくことであ る。以下がそれを実現するプログラムになる。
pt minus 1 pt
void bhnode::assign_root(vector root_pos, real length, particle * p, int np)
{
pos = root_pos;
l = length;
pfirst = p;
nparticle = np;
for(int i=0;i<np-1;i++){
p->next = p+1;
p++;
}
p->next = NULL;
}
これはメンバ関数なので、bhnode クラス変数のメンバはメンバ名だけ
で使える。中心座標、辺の長さと、粒子の配列、粒子数が引数で渡ってくるの
で、中心座標、辺の長さをメンバ変数に代入し、 pfirst から始まるリ
ストを配列の要素を順にたどっていくように作っておく。
次に実際に木構造を作っていく。これは、以下のようになる
pt minus 1 pt
void bhnode::create_tree_recursive(bhnode * & heap_top, int & heap_remainder)
{
if (heap_remainder <= 0){
cerr << "create_tree: no more free node... exit\n";
exit(1);
}
particle * psub[8];
int nsub[8];
for(int i = 0;i<8;i++){
psub[i] = NULL;
nsub[i] = 0;
}
particle * p = pfirst;
for(int i=0 ; i<nparticle; i++, p = p->next){
int subindex = 0;
for(int k=0;k<3;k++){
subindex <<= 1;
if (p->pos[k] > cpos[k]){
subindex += 1;
}
}
nsub[subindex] += 1;
p->subnext = psub[subindex];
psub[subindex] = p;
}
for(int i=0; i<8;i++) if(nsub[i] >0){
child[i] = heap_top;
heap_top ++;
heap_remainder -- ;
child[i]->pfirst = psub[i];
child[i]->cpos = cpos + vector( ((i&4)*0.5-1)*l/4,
((i&2) -1)*l/4,
((i&1)*2 -1)*l/4);
child[i]->l = l*0.5;
child[i]->nparticle = nsub[i];
for(p=psub[i]; p->subnext != NULL; p = p->subnext){
p->next = p->subnext;
}
if (nsub[i] > 1){
child[i]->create_tree_recursive(heap_top, heap_remainder);
}else{
child[i]->pos = child[i]->pfirst->pos;
child[i]->mass = child[i]->pfirst->mass;
}
}
}
ここで、渡ってきている引数は、実際には bhnode 構造体の配列の先頭アド
レスと、その配列の要素数、というか正確にはその配列のどこか途中のアドレ
スと残りの要素数である。この配列は、木構造全部を作るのに十分な程度の大
きさをとっておく。木に新しいノードを付け加えるためには配列の先頭から順
番に要素を割り当てていく。
一般には、木構造のような動的な構造を作るには、 new 演算子を使っ て必要に応じて新しいメモリ領域を割り当てる。が、ここでは木構造は時間積分 のステップ毎に新しく作るので、 new をその度に使ってはメモリが足 りなくなるとかややこしい問題がある。そこで、自前でメモリを配列に確保してお いてそこにあるものを使うことにする。足りなくなったらここではエラーを出 して終る。
まず、最初にすることは、自分の下にある全粒子をたどっていって、それがど の子ノードにいくかを判定し、結果に応じて対応するリストに挿入（先頭に追 加）することである。ここで、 next は自分のリストなので、これを書 き変えるとその後がたどれなくなる。このため、子ノード用のリストには subnext を使うことにする。判定は、各座標方向に対して、中心より上 か下かを見るだけである。
全部の粒子の振り分けがすんだら、次は1つ以上の粒子が入っている子ノード にちゃんとクラス変数を割り当て、その中身（位置と大きさ）を設定するのと 同時にリストのほうも subnextを next にコピーしなおしておく。
で、粒子の数が 2個以上なら、今度はこの子ノードをさらに分割する。これは、 ここではこの関数自身を再帰的に呼び出すことで実現されている。粒子が1つ だけなら、その質量、位置を木構造データのほうの変数にコピーしておく。
と、今日はとりあえずこれくらいで。来週はこれを使った重力の計算をしてみ る。