| 1 | /****************************************************************************/ |
|---|
| 2 | /* TREELOAD.C: routines to create tree. Public routines: maketree(). */ |
|---|
| 3 | /* Copyright (c) 1999 by Joshua E. Barnes, Tokyo, JAPAN. */ |
|---|
| 4 | /****************************************************************************/ |
|---|
| 5 | |
|---|
| 6 | #include "stdinc.h" |
|---|
| 7 | #include "mathfns.h" |
|---|
| 8 | #include "vectmath.h" |
|---|
| 9 | #include "treedefs.h" |
|---|
| 10 | |
|---|
| 11 | /* |
|---|
| 12 | * Local routines and variables to perform tree construction. |
|---|
| 13 | */ |
|---|
| 14 | |
|---|
| 15 | local void newtree(void); /* flush existing tree */ |
|---|
| 16 | local cellptr makecell(void); /* create an empty cell */ |
|---|
| 17 | local void expandbox(bodyptr, int); /* set size of root cell */ |
|---|
| 18 | local void loadbody(bodyptr); /* load body into tree */ |
|---|
| 19 | local int subindex(bodyptr, cellptr); /* compute subcell index */ |
|---|
| 20 | local void hackcofm(cellptr, real, int); /* find centers of mass */ |
|---|
| 21 | local void setrcrit(cellptr, vector, real); /* set cell's crit. radius */ |
|---|
| 22 | local void threadtree(nodeptr, nodeptr); /* set next and more links */ |
|---|
| 23 | local void hackquad(cellptr); /* compute quad moments */ |
|---|
| 24 | |
|---|
| 25 | local bool bh86, sw94; /* use alternate criteria */ |
|---|
| 26 | local nodeptr freecell = NULL; /* list of free cells */ |
|---|
| 27 | |
|---|
| 28 | #define MAXLEVEL 32 /* max height of tree */ |
|---|
| 29 | |
|---|
| 30 | local int cellhist[MAXLEVEL]; /* count cells by level */ |
|---|
| 31 | local int subnhist[MAXLEVEL]; /* count subnodes by level */ |
|---|
| 32 | |
|---|
| 33 | /* |
|---|
| 34 | * MAKETREE: initialize tree structure for hierarchical force calculation. |
|---|
| 35 | */ |
|---|
| 36 | |
|---|
| 37 | void maketree(bodyptr btab, int nbody) |
|---|
| 38 | { |
|---|
| 39 | double cpustart; |
|---|
| 40 | bodyptr p; |
|---|
| 41 | int i; |
|---|
| 42 | |
|---|
| 43 | cpustart = cputime(); /* record time at start */ |
|---|
| 44 | newtree(); /* flush existing tree, etc */ |
|---|
| 45 | root = makecell(); /* allocate the root cell */ |
|---|
| 46 | CLRV(Pos(root)); /* initialize the midpoint */ |
|---|
| 47 | expandbox(btab, nbody); /* and expand cell to fit */ |
|---|
| 48 | for (p = btab; p < btab+nbody; p++) /* loop over all bodies */ |
|---|
| 49 | loadbody(p); /* insert each into tree */ |
|---|
| 50 | bh86 = scanopt(options, "bh86"); /* set flags for alternate */ |
|---|
| 51 | sw94 = scanopt(options, "sw94"); /* ...cell opening criteria */ |
|---|
| 52 | if (bh86 && sw94) /* can't have both at once */ |
|---|
| 53 | error("maketree: incompatible options bh86 and sw94\n"); |
|---|
| 54 | tdepth = 0; /* init count of levels */ |
|---|
| 55 | for (i = 0; i < MAXLEVEL; i++) /* and init tree histograms */ |
|---|
| 56 | cellhist[i] = subnhist[i] = 0; |
|---|
| 57 | hackcofm(root, rsize, 0); /* find c-of-m coords, etc */ |
|---|
| 58 | threadtree((nodeptr) root, NULL); /* add next and more links */ |
|---|
| 59 | if (usequad) /* if including quad terms */ |
|---|
| 60 | hackquad(root); /* find quadrupole moments */ |
|---|
| 61 | cputree = cputime() - cpustart; /* store elapsed CPU time */ |
|---|
| 62 | } |
|---|
| 63 | |
|---|
| 64 | /* |
|---|
| 65 | * NEWTREE: reclaim cells in tree, prepare to build new one. |
|---|
| 66 | */ |
|---|
| 67 | |
|---|
| 68 | local void newtree(void) |
|---|
| 69 | { |
|---|
| 70 | static bool firstcall = TRUE; |
|---|
| 71 | nodeptr p; |
|---|
| 72 | |
|---|
| 73 | if (! firstcall) { /* if cells to reclaim */ |
|---|
| 74 | p = (nodeptr) root; /* start with the root */ |
|---|
| 75 | while (p != NULL) /* loop scanning tree */ |
|---|
| 76 | if (Type(p) == CELL) { /* if we found a cell to */ |
|---|
| 77 | Next(p) = freecell; /* then save existing list */ |
|---|
| 78 | freecell = p; /* and add it to the front */ |
|---|
| 79 | p = More(p); /* then scan down tree */ |
|---|
| 80 | } else /* else, skip over bodies */ |
|---|
| 81 | p = Next(p); /* by going on to the next */ |
|---|
| 82 | } else /* else nothing to reclaim */ |
|---|
| 83 | firstcall = FALSE; /* so just note it */ |
|---|
| 84 | root = NULL; /* flush existing tree */ |
|---|
| 85 | ncell = 0; /* reset cell count */ |
|---|
| 86 | } |
|---|
| 87 | |
|---|
| 88 | /* |
|---|
| 89 | * MAKECELL: return pointer to free cell. |
|---|
| 90 | */ |
|---|
| 91 | |
|---|
| 92 | local cellptr makecell(void) |
|---|
| 93 | { |
|---|
| 94 | cellptr c; |
|---|
| 95 | int i; |
|---|
| 96 | |
|---|
| 97 | if (freecell == NULL) /* if no free cells left */ |
|---|
| 98 | c = (cellptr) allocate(sizeof(cell)); /* then allocate a new one */ |
|---|
| 99 | else { /* else use existing cell */ |
|---|
| 100 | c = (cellptr) freecell; /* take the one in front */ |
|---|
| 101 | freecell = Next(c); /* and go on to next one */ |
|---|
| 102 | } |
|---|
| 103 | Type(c) = CELL; /* initialize node type */ |
|---|
| 104 | Update(c) = FALSE; /* and force update flag */ |
|---|
| 105 | for (i = 0; i < NSUB; i++) /* loop over subcells */ |
|---|
| 106 | Subp(c)[i] = NULL; /* and empty each one */ |
|---|
| 107 | ncell++; /* count one more cell */ |
|---|
| 108 | return (c); /* return pointer to cell */ |
|---|
| 109 | } |
|---|
| 110 | |
|---|
| 111 | /* |
|---|
| 112 | * EXPANDBOX: find range of coordinate values (with respect to root) |
|---|
| 113 | * and expand root cell to fit. The size is doubled at each step to |
|---|
| 114 | * take advantage of exact representation of powers of two. |
|---|
| 115 | */ |
|---|
| 116 | |
|---|
| 117 | local void expandbox(bodyptr btab, int nbody) |
|---|
| 118 | { |
|---|
| 119 | real dmax, d; |
|---|
| 120 | bodyptr p; |
|---|
| 121 | int k; |
|---|
| 122 | |
|---|
| 123 | dmax = 0.0; /* keep track of max value */ |
|---|
| 124 | for (p = btab; p < btab+nbody; p++) /* loop over all bodies */ |
|---|
| 125 | for (k = 0; k < NDIM; k++) { /* and over all dimensions */ |
|---|
| 126 | d = rabs(Pos(p)[k] - Pos(root)[k]); /* find distance to midpnt */ |
|---|
| 127 | if (d > dmax) /* if bigger than old one */ |
|---|
| 128 | dmax = d; /* store new max value */ |
|---|
| 129 | } |
|---|
| 130 | while (rsize < 2 * dmax) /* loop until value fits */ |
|---|
| 131 | rsize = 2 * rsize; /* doubling box each time */ |
|---|
| 132 | } |
|---|
| 133 | |
|---|
| 134 | /* |
|---|
| 135 | * LOADBODY: descend tree and insert body p in appropriate place. |
|---|
| 136 | */ |
|---|
| 137 | |
|---|
| 138 | local void loadbody(bodyptr p) |
|---|
| 139 | { |
|---|
| 140 | cellptr q, c; |
|---|
| 141 | int qind, k; |
|---|
| 142 | real qsize, dist2; |
|---|
| 143 | vector distv; |
|---|
| 144 | |
|---|
| 145 | q = root; /* start with tree root */ |
|---|
| 146 | qind = subindex(p, q); /* get index of subcell */ |
|---|
| 147 | qsize = rsize; /* keep track of cell size */ |
|---|
| 148 | while (Subp(q)[qind] != NULL) { /* loop descending tree */ |
|---|
| 149 | if (Type(Subp(q)[qind]) == BODY) { /* if another body reached */ |
|---|
| 150 | DOTPSUBV(dist2, distv, Pos(p), Pos(Subp(q)[qind])); |
|---|
| 151 | if (dist2 == 0.0) /* check positions differ */ |
|---|
| 152 | error("loadbody: two bodies have same position\n"); |
|---|
| 153 | c = makecell(); /* then allocate new cell */ |
|---|
| 154 | for (k = 0; k < NDIM; k++) /* and initialize midpoint */ |
|---|
| 155 | Pos(c)[k] = Pos(q)[k] + /* offset from parent */ |
|---|
| 156 | (Pos(p)[k] < Pos(q)[k] ? - qsize : qsize) / 4; |
|---|
| 157 | Subp(c)[subindex((bodyptr) Subp(q)[qind], c)] = Subp(q)[qind]; |
|---|
| 158 | /* put body in cell */ |
|---|
| 159 | Subp(q)[qind] = (nodeptr) c; /* link cell in tree */ |
|---|
| 160 | } |
|---|
| 161 | q = (cellptr) Subp(q)[qind]; /* advance to next level */ |
|---|
| 162 | qind = subindex(p, q); /* get index to examine */ |
|---|
| 163 | qsize = qsize / 2; /* shrink current cell */ |
|---|
| 164 | } |
|---|
| 165 | Subp(q)[qind] = (nodeptr) p; /* found place, store p */ |
|---|
| 166 | } |
|---|
| 167 | |
|---|
| 168 | /* |
|---|
| 169 | * SUBINDEX: compute subcell index for body p in cell q. |
|---|
| 170 | */ |
|---|
| 171 | |
|---|
| 172 | local int subindex(bodyptr p, cellptr q) |
|---|
| 173 | { |
|---|
| 174 | int ind, k; |
|---|
| 175 | |
|---|
| 176 | ind = 0; /* accumulate subcell index */ |
|---|
| 177 | for (k = 0; k < NDIM; k++) /* loop over dimensions */ |
|---|
| 178 | if (Pos(q)[k] <= Pos(p)[k]) /* if beyond midpoint */ |
|---|
| 179 | ind += NSUB >> (k + 1); /* then skip over subcells */ |
|---|
| 180 | return (ind); |
|---|
| 181 | } |
|---|
| 182 | |
|---|
| 183 | /* |
|---|
| 184 | * HACKCOFM: descend tree finding center-of-mass coordinates; |
|---|
| 185 | * also sets critical cell radii, if appropriate. |
|---|
| 186 | */ |
|---|
| 187 | |
|---|
| 188 | local void hackcofm(cellptr p, real psize, int lev) |
|---|
| 189 | { |
|---|
| 190 | vector cmpos, tmpv; |
|---|
| 191 | int i, k; |
|---|
| 192 | nodeptr q; |
|---|
| 193 | real dpq; |
|---|
| 194 | |
|---|
| 195 | tdepth = MAX(tdepth, lev); /* remember maximum level */ |
|---|
| 196 | cellhist[lev]++; /* count cells by level */ |
|---|
| 197 | Mass(p) = 0.0; /* init cell's total mass */ |
|---|
| 198 | CLRV(cmpos); /* and center of mass pos */ |
|---|
| 199 | for (i = 0; i < NSUB; i++) /* loop over the subnodes */ |
|---|
| 200 | if ((q = Subp(p)[i]) != NULL) { /* skipping the NULLs */ |
|---|
| 201 | subnhist[lev]++; /* count existing subnodes */ |
|---|
| 202 | if (Type(q) == CELL) /* and if node is a cell */ |
|---|
| 203 | hackcofm((cellptr) q, psize/2, lev+1); |
|---|
| 204 | /* then do the same for it */ |
|---|
| 205 | Update(p) |= Update(q); /* propagate update request */ |
|---|
| 206 | Mass(p) += Mass(q); /* accumulate total mass */ |
|---|
| 207 | MULVS(tmpv, Pos(q), Mass(q)); /* weight position by mass */ |
|---|
| 208 | ADDV(cmpos, cmpos, tmpv); /* and sum c-of-m position */ |
|---|
| 209 | } |
|---|
| 210 | if (Mass(p) > 0.0) { /* usually, cell has mass */ |
|---|
| 211 | DIVVS(cmpos, cmpos, Mass(p)); /* so find c-of-m position */ |
|---|
| 212 | } else { /* but if no mass inside */ |
|---|
| 213 | SETV(cmpos, Pos(p)); /* use geo. center for now */ |
|---|
| 214 | } |
|---|
| 215 | for (k = 0; k < NDIM; k++) /* check c-of-m of cell */ |
|---|
| 216 | if (cmpos[k] < Pos(p)[k] - psize/2 || /* if actually outside cell */ |
|---|
| 217 | Pos(p)[k] + psize/2 <= cmpos[k]) |
|---|
| 218 | error("hackcofm: tree structure error\n"); |
|---|
| 219 | #if !defined(QUICKSCAN) |
|---|
| 220 | setrcrit(p, cmpos, psize); /* set critical radius */ |
|---|
| 221 | #endif |
|---|
| 222 | SETV(Pos(p), cmpos); /* and center-of-mass pos */ |
|---|
| 223 | } |
|---|
| 224 | |
|---|
| 225 | #if !defined(QUICKSCAN) |
|---|
| 226 | |
|---|
| 227 | /* |
|---|
| 228 | * SETRCRIT: assign critical radius for cell p, using center-of-mass |
|---|
| 229 | * position cmpos and cell size psize. |
|---|
| 230 | */ |
|---|
| 231 | |
|---|
| 232 | local void setrcrit(cellptr p, vector cmpos, real psize) |
|---|
| 233 | { |
|---|
| 234 | real bmax2, d; |
|---|
| 235 | int k; |
|---|
| 236 | |
|---|
| 237 | if (theta == 0.0) /* if exact calculation */ |
|---|
| 238 | Rcrit2(p) = rsqr(2 * rsize); /* then always open cells */ |
|---|
| 239 | else if (sw94) { /* if using S&W's criterion */ |
|---|
| 240 | bmax2 = 0.0; /* compute max distance^2 */ |
|---|
| 241 | for (k = 0; k < NDIM; k++) { /* loop over dimensions */ |
|---|
| 242 | d = cmpos[k] - Pos(p)[k] + psize/2; /* get dist from corner */ |
|---|
| 243 | bmax2 += rsqr(MAX(d, psize - d)); /* and sum max distance^2 */ |
|---|
| 244 | } |
|---|
| 245 | Rcrit2(p) = bmax2 / rsqr(theta); /* use max dist from cm */ |
|---|
| 246 | } else if (bh86) /* if using old criterion */ |
|---|
| 247 | Rcrit2(p) = rsqr(psize / theta); /* then use size of cell */ |
|---|
| 248 | else { /* else use new criterion */ |
|---|
| 249 | DISTV(d, cmpos, Pos(p)); /* find offset from center */ |
|---|
| 250 | Rcrit2(p) = rsqr(psize / theta + d); /* use size plus offset */ |
|---|
| 251 | } |
|---|
| 252 | } |
|---|
| 253 | |
|---|
| 254 | #endif |
|---|
| 255 | |
|---|
| 256 | /* |
|---|
| 257 | * THREADTREE: do a recursive treewalk starting from node p, |
|---|
| 258 | * with next stop n, installing Next and More links. |
|---|
| 259 | */ |
|---|
| 260 | |
|---|
| 261 | local void threadtree(nodeptr p, nodeptr n) |
|---|
| 262 | { |
|---|
| 263 | int ndesc, i; |
|---|
| 264 | nodeptr desc[NSUB+1]; |
|---|
| 265 | |
|---|
| 266 | Next(p) = n; /* set link to next node */ |
|---|
| 267 | if (Type(p) == CELL) { /* if descendents to thread */ |
|---|
| 268 | ndesc = 0; /* start counting them */ |
|---|
| 269 | for (i = 0; i < NSUB; i++) /* loop over all subcells */ |
|---|
| 270 | if (Subp(p)[i] != NULL) /* if this one is occupied */ |
|---|
| 271 | desc[ndesc++] = Subp(p)[i]; /* then store it in table */ |
|---|
| 272 | More(p) = desc[0]; /* set link to 1st one */ |
|---|
| 273 | desc[ndesc] = n; /* thread last one to next */ |
|---|
| 274 | for (i = 0; i < ndesc; i++) /* loop over descendents */ |
|---|
| 275 | threadtree(desc[i], desc[i+1]); /* and thread them together */ |
|---|
| 276 | } |
|---|
| 277 | } |
|---|
| 278 | |
|---|
| 279 | /* |
|---|
| 280 | * HACKQUAD: descend tree, evaluating quadrupole moments. Note that this |
|---|
| 281 | * routine is coded so that the Subp() and Quad() components of a cell can |
|---|
| 282 | * share the same memory locations. |
|---|
| 283 | */ |
|---|
| 284 | |
|---|
| 285 | local void hackquad(cellptr p) |
|---|
| 286 | { |
|---|
| 287 | int ndesc, i; |
|---|
| 288 | nodeptr desc[NSUB], q; |
|---|
| 289 | vector dr; |
|---|
| 290 | real drsq; |
|---|
| 291 | matrix drdr, Idrsq, tmpm; |
|---|
| 292 | |
|---|
| 293 | ndesc = 0; /* count occupied subnodes */ |
|---|
| 294 | for (i = 0; i < NSUB; i++) /* loop over all subnodes */ |
|---|
| 295 | if (Subp(p)[i] != NULL) /* if this one's occupied */ |
|---|
| 296 | desc[ndesc++] = Subp(p)[i]; /* copy it to safety */ |
|---|
| 297 | CLRM(Quad(p)); /* init quadrupole moment */ |
|---|
| 298 | for (i = 0; i < ndesc; i++) { /* loop over real subnodes */ |
|---|
| 299 | q = desc[i]; /* access ech one in turn */ |
|---|
| 300 | if (Type(q) == CELL) /* if it's also a cell */ |
|---|
| 301 | hackquad((cellptr) q); /* then process it first */ |
|---|
| 302 | SUBV(dr, Pos(q), Pos(p)); /* find displacement vect. */ |
|---|
| 303 | OUTVP(drdr, dr, dr); /* form outer prod. of dr */ |
|---|
| 304 | DOTVP(drsq, dr, dr); /* and dot prod. dr * dr */ |
|---|
| 305 | SETMI(Idrsq); /* init unit matrix */ |
|---|
| 306 | MULMS(Idrsq, Idrsq, drsq); /* and scale by dr * dr */ |
|---|
| 307 | MULMS(tmpm, drdr, 3.0); /* scale drdr by 3 */ |
|---|
| 308 | SUBM(tmpm, tmpm, Idrsq); /* now form quad. moment */ |
|---|
| 309 | MULMS(tmpm, tmpm, Mass(q)); /* from cm of subnode */ |
|---|
| 310 | if (Type(q) == CELL) /* if subnode is cell */ |
|---|
| 311 | ADDM(tmpm, tmpm, Quad(q)); /* then include its moment */ |
|---|
| 312 | ADDM(Quad(p), Quad(p), tmpm); /* increment moment of cell */ |
|---|
| 313 | } |
|---|
| 314 | } |
|---|
![(please configure the [header_logo] section in trac.ini)](/trac/PP2009/chrome/site/your_project_logo.png){kind=logo}
