1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
|
/****************************************************************************/
/* LOAD.C: routines to create tree. Public routines: maketree(). */
/* */
/* Copyright (c) 1993 by Joshua E. Barnes, Honolulu, HI. */
/* It's free because it's yours. */
/****************************************************************************/
#include "defs.h"
local void newtree(void); /* flush existing tree */
local cellptr makecell(void); /* create an empty cell */
local void expandbox(bodyptr, int); /* set size of root cell */
local void loadbody(bodyptr); /* load body into tree */
local int subindex(bodyptr, cellptr); /* compute subcell index */
local void hackcofm(cellptr, real); /* find centers of mass */
local void setrcrit(cellptr, vector, real); /* set cell's crit. radius */
local void threadtree(nodeptr, nodeptr); /* set next and more links */
local void hackquad(cellptr); /* compute quad moments */
/*
* MAKETREE: initialize tree structure for hierarchical force calculation
* from body array btab, which contains nbody bodies.
*/
local bool bh86, sw93; /* use alternate criteria */
void maketree(bodyptr btab, int nbody)
{
bodyptr p;
newtree(); /* flush existing tree, etc */
root = makecell(); /* allocate the root cell */
CLRV(Pos(root)); /* initialize the midpoint */
expandbox(btab, nbody); /* and expand cell to fit */
maxlevel = 0; /* init count of levels */
for (p = btab; p < btab+nbody; p++) /* loop over bodies... */
if (Mass(p) != 0.0) /* exclude test particles */
loadbody(p); /* and insert into tree */
bh86 = scanopt(options, "bh86"); /* set flags for alternate */
sw93 = scanopt(options, "sw93"); /* ...cell opening criteria */
if (bh86 && sw93) /* can't have both at once */
error("maketree: options bh86 and sw93 are incompatible\n");
hackcofm(root, rsize); /* find c-of-m coordinates */
threadtree((nodeptr) root, NULL); /* add Next and More links */
if (usequad) /* including quad moments? */
hackquad(root); /* assign Quad moments */
}
/*
* NEWTREE: reclaim cells in tree, prepare to build new one.
*/
local nodeptr freecell = NULL; /* list of free cells */
local void newtree(void)
{
permanent bool firstcall = TRUE;
nodeptr p;
if (! firstcall) { /* tree data to reclaim? */
p = (nodeptr) root; /* start with the root */
while (p != NULL) /* loop scanning tree */
if (Type(p) == CELL) { /* found cell to free? */
Next(p) = freecell; /* link to front of */
freecell = p; /* ...existing list */
p = MoreN(p); /* scan down tree */
} else /* skip over bodies */
p = Next(p); /* go on to next */
} else /* first time through */
firstcall = FALSE; /* so just note it */
root = NULL; /* flush existing tree */
cellused = 0; /* reset cell count */
}
/*
* MAKECELL: return pointer to free cell.
*/
local cellptr makecell(void)
{
cellptr c;
int i;
if (freecell == NULL) { /* no free cells left? */
c = (cellptr) allocate(sizeof(cell)); /* allocate a new one */
#ifndef NO_PERF_CHANGES
Type(c) = CELL; /* initialize cell type */
#endif
} else { /* use existing free cell */
Type(freecell) = CELL; /* initialize cell type */
c = node2cell (freecell); /* take one on front */
freecell = Next(c); /* go on to next one */
}
#ifdef NO_PERF_CHANGES
Type(c) = CELL; /* initialize cell type */
#endif
for (i = 0; i < NSUB; i++) /* loop over subcells */
Subp(c)[i] = NULL; /* and empty each one */
cellused++; /* count one more cell */
return (c);
}
/*
* EXPANDBOX: find range of coordinate values (with respect to root)
* and expand root cell to fit. The size is doubled at each step to
* take advantage of exact representation of powers of two.
*/
local void expandbox(bodyptr btab, int nbody)
{
real xyzmax;
bodyptr p;
int k;
xyzmax = 0.0;
for (p = btab; p < btab+nbody; p++)
for (k = 0; k < NDIM; k++)
xyzmax = MAX(xyzmax, rabs(Pos(p)[k] - Pos(root)[k]));
while (rsize < 2 * xyzmax)
rsize = 2 * rsize;
}
/*
* LOADBODY: descend tree and insert body p in appropriate place.
*/
local void loadbody(bodyptr p)
{
cellptr q, c;
int qind, lev, k;
real qsize;
q = root; /* start with tree root */
qind = subindex(p, q); /* get index of subcell */
qsize = rsize; /* keep track of cell size */
lev = 0; /* count levels descended */
while (Subp(q)[qind] != NULL) { /* loop descending tree */
if (Type(Subp(q)[qind]) == BODY) { /* reached a "leaf"? */
c = makecell(); /* allocate new cell */
for (k = 0; k < NDIM; k++) /* initialize midpoint */
Pos(c)[k] = Pos(q)[k] + /* offset from parent */
(Pos(p)[k] < Pos(q)[k] ? - qsize : qsize) / 4;
Subp(c)[subindex(node2body(Subp(q)[qind]), c)] = Subp(q)[qind];
/* put body in cell */
Subp(q)[qind] = (nodeptr) c; /* link cell in tree */
}
q = node2cell(Subp(q)[qind]); /* advance to next level */
qind = subindex(p, q); /* get index to examine */
qsize = qsize / 2; /* shrink current cell */
lev++; /* count another level */
}
Subp(q)[qind] = (nodeptr) p; /* found place, store p */
maxlevel = MAX(maxlevel, lev); /* remember maximum level */
}
/*
* SUBINDEX: compute subcell index for body p in cell q.
*/
local int subindex(bodyptr p, cellptr q)
{
int ind, k;
ind = 0; /* accumulate subcell index */
for (k = 0; k < NDIM; k++) /* loop over dimensions */
if (Pos(q)[k] <= Pos(p)[k]) /* if beyond midpoint */
ind += NSUB >> (k + 1); /* skip over subcells */
return (ind);
}
/*
* HACKCOFM: descend tree finding center-of-mass coordinates and
* setting critical cell radii.
*/
local void hackcofm(cellptr p, real psize)
{
vector cmpos, tmpv;
int i, k;
nodeptr q;
Mass(p) = 0.0; /* init total mass... */
CLRV(cmpos); /* and center of mass */
for (i = 0; i < NSUB; i++) /* loop over subnodes */
if ((q = Subp(p)[i]) != NULL) { /* does subnode exist? */
if (Type(q) == CELL) /* and is it a cell? */
hackcofm(node2cell(q), psize/2); /* find subcell cm */
Mass(p) += Mass(q); /* sum total mass */
MULVS(tmpv, Pos(q), Mass(q)); /* weight pos by mass */
ADDV(cmpos, cmpos, tmpv); /* sum c-of-m position */
}
DIVVS(cmpos, cmpos, Mass(p)); /* rescale cms position */
for (k = 0; k < NDIM; k++) /* check tree structure... */
if (cmpos[k] < Pos(p)[k] - psize/2 || /* if out of bounds */
Pos(p)[k] + psize/2 <= cmpos[k]) /* in either direction */
error("hackcofm: tree structure error\n");
setrcrit(p, cmpos, psize); /* set critical radius */
SETV(Pos(p), cmpos); /* and center-of-mass pos */
}
/*
* SETRCRIT: assign critical radius for cell p, using center-of-mass
* position cmpos and cell size psize.
*/
local void setrcrit(cellptr p, vector cmpos, real psize)
{
real rc, bmax2, dmin;
int k;
if (theta == 0.0) /* exact force calculation? */
rc = 2 * rsize; /* always open cells */
else if (bh86) /* use old BH criterion? */
rc = psize / theta; /* using size of cell */
else if (sw93) { /* use S&W's criterion? */
bmax2 = 0.0; /* compute max distance^2 */
for (k = 0; k < NDIM; k++) { /* loop over dimensions */
dmin = cmpos[k] - (Pos(p)[k] - psize/2);
/* dist from 1st corner */
bmax2 += rsqr(MAX(dmin, psize - dmin));
/* sum max distance^2 */
}
rc = rsqrt(bmax2) / theta; /* using max dist from cm */
} else { /* use new criterion? */
rc = psize / theta + distv(cmpos, Pos(p));
/* use size plus offset */
}
Rcrit2(p) = rsqr(rc); /* store square of radius */
}
/*
* THREADTREE: do a recursive treewalk starting from node p,
* with next stop n, installing Next and More links.
*/
local void threadtree(nodeptr p, nodeptr n)
{
int ndesc, i;
nodeptr desc[NSUB+1];
Next(p) = n; /* link to next node */
if (Type(p) == CELL) { /* any children to thread? */
ndesc = 0; /* count extant children */
for (i = 0; i < NSUB; i++) /* loop over subnodes */
if (SubpN(p)[i] != NULL) /* found a live one? */
desc[ndesc++] = SubpN(p)[i]; /* store in table */
MoreN(p) = desc[0]; /* link to first child */
desc[ndesc] = n; /* end table with next */
for (i = 0; i < ndesc; i++) /* loop over children */
threadtree(desc[i], desc[i+1]); /* thread each w/ next */
}
}
/*
* HACKQUAD: descend tree, evaluating quadrupole moments. Note that this
* routine is coded so that the Subp() and Quad() components of a cell can
* share the same memory locations.
*/
local void hackquad(cellptr p)
{
int i;
nodeptr psub[NSUB], q;
vector dr;
real drsq;
matrix drdr, Idrsq, tmpm;
for (i = 0; i < NSUB; i++) /* loop over subnodes */
psub[i] = Subp(p)[i]; /* copy each to safety */
CLRM(Quad(p)); /* init quadrupole moment */
for (i = 0; i < NSUB; i++) /* loop over subnodes */
if ((q = psub[i]) != NULL) { /* does subnode exist? */
if (Type(q) == CELL) /* and is it a call? */
hackquad(node2cell(q)); /* process it first */
SUBV(dr, Pos(q), Pos(p)); /* displacement vect. */
OUTVP(drdr, dr, dr); /* outer prod. of dr */
DOTVP(drsq, dr, dr); /* dot prod. dr * dr */
SETMI(Idrsq); /* init unit matrix */
MULMS(Idrsq, Idrsq, drsq); /* scale by dr * dr */
MULMS(tmpm, drdr, 3.0); /* scale drdr by 3 */
SUBM(tmpm, tmpm, Idrsq); /* form quad. moment */
MULMS(tmpm, tmpm, Mass(q)); /* from cm of subnode */
if (Type(q) == CELL) /* if subnode is cell */
ADDM(tmpm, tmpm, QuadN(q)); /* add its moment */
ADDM(Quad(p), Quad(p), tmpm); /* add to qm of cell */
}
}
|
