// // BHtree.C #define PR(x) cerr << #x << " = " << x << " " #define PRC(x) cerr << #x << " = " << x << ", " #define PRL(x) cerr << #x << " = " << x << "\n" #include #include #include #define real double #include "BHtree.h" void bhnode::assign_root(myvector root_pos, real length, particle * p, int np) { pos = root_pos; l = length; pfirst = p; nparticle = np; for(int i=0;inext = p+1; p++; } p->next = NULL; } int childindex(myvector pos, myvector cpos) { int subindex=0; for(int k=0;k<3;k++){ subindex <<= 1; if (pos[k] > cpos[k]){ subindex += 1; } } return subindex; } void bhnode::assign_child(int subindex, bhnode * & heap_top, int & heap_remainder) { if (heap_remainder <= 0){ cerr << "create_tree: no more free node... exit\n"; exit(1); } if (child[subindex]==NULL){ child[subindex] = heap_top; heap_top ++; heap_remainder -- ; child[subindex]->cpos = cpos + myvector( ((subindex&4)*0.5-1)*l/4, ((subindex&2) -1)*l/4, ((subindex&1)*2 -1)*l/4); child[subindex]->l = l*0.5; child[subindex]->nparticle=0; } } void bhnode::create_tree_recursive(bhnode * & heap_top, int & heap_remainder) { for(int i = 0;i<8;i++)child[i] = NULL; particle * p = pfirst; for(int i=0 ; inext; int subindex = childindex(p->pos,cpos); assign_child(subindex, heap_top, heap_remainder); child[subindex]->nparticle ++; p->next = child[subindex]->pfirst; child[subindex]->pfirst = p; p=pnext; } for(int i=0; i<8;i++) if(child[i]!=NULL){ if (child[i]->nparticle > 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; } } } void spc(int indent) { for(int i=0;ipos); p= p->next; } }else{ cerr << " IS _not_ LEAF ";PRL(nparticle); for(i=0;i<8;i++){ if (child[i] != NULL){ child[i]->dump(indent + 2); } } } } void bhnode::dumptree(int indent) { int i; spc(indent); cerr << "node center pos " << cpos ; if (nparticle == 1){ cerr << " IS LEAF" ;PRL(nparticle); particle * p = pfirst; for(i = 0; i < nparticle; i++){ for(int j=0;jpos); p= p->next; } }else{ cerr << " IS _not_ LEAF ";PRL(nparticle); for(i=0;i<8;i++){ if (child[i] != NULL){ child[i]->dumptree(indent + 2); } } } } int inbox(myvector & cpos, // center of the box myvector & pos, // position of the particle real l) // length of one side of the box { for(int i = 0; i< ndim; i++){ if (fabs(pos[i]-cpos[i]) > l*0.5) return 1; } return 0; } int bhnode::sanity_check() { int i; int iret = 0; if (nparticle == 1 ){ // this is the lowest level node. Things to check: // particle is in the cell particle * p = pfirst; if(inbox(pos,p->pos,l)){ cerr << "Error, particle out of box ... \n"; dump(); return 1; } }else{ // This is the non-leaf node. Check the position and side // length of the child cells and then check recursively.. for(i=0;i<8;i++){ if (child[i] != NULL){ int err = 0; err = child[i]->sanity_check(); if (l*0.5 != child[i]->l) err += 2; myvector relpos = cpos-child[i]->cpos; for (int k = 0 ; k> n ; pp = new particle[n]; double rsize = 1.0; create_uniform_sphere(pp, n, 0 , rsize); for(int i =0;i assign_root(myvector(0.0), rsize*2, pp, n); int heap_remainder = nnodes-1; bhnode * btmp = bn + 1; bn->create_tree_recursive(btmp, heap_remainder); PRL(bn->sanity_check()); bn->set_cm_quantities(); bn->dump(); double eps2 = 0.01; calculate_uncorrected_gravity_direct(pp, n, eps2); cerr << "Direct force \n"; particle * p = pp; for(int i = 0; ipos);PR(p->phi); PRL(p->acc); p++; } cerr << "Tree force \n"; clear_acc_and_phi(pp, n); p = pp; for(int i = 0; ipos);PR(p->phi); PRL(p->acc); p++; } return 0; } #endif #ifdef TREETEST int main() { particle * pp; int n; cerr << "Enter n:"; cin >> n ; pp = new particle[n]; double rsize = 1.0; create_uniform_sphere(pp, n, 0 , rsize); for(int i =0;i assign_root(myvector(0.0), rsize*2, pp, n); int heap_remainder = nnodes-1; bhnode * btmp = bn + 1; bn->create_tree_recursive(btmp, heap_remainder); PRL(bn->sanity_check()); bn->dumptree(); return 0; } #endif