Timothy  0.9
Tissue Modelling Framework
 All Data Structures Files Functions Variables Typedefs Macros
tree.c
Go to the documentation of this file.
1 /* **************************************************************************
2  * This file is part of Timothy
3  *
4  * Copyright (c) 2014/15 Maciej Cytowski
5  * Copyright (c) 2014/15 ICM, University of Warsaw, Poland
6  *
7  * This program is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 2 of the License, or
10  * (at your option) any later version.
11  *
12  * This program is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with this program; if not, write to the Free Software
19  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
20  *
21  * *************************************************************************/
22 
23 #include <stdio.h>
24 #include <stdlib.h>
25 #include <math.h>
26 
27 #include "global.h"
28 #include "inline.h"
29 
34 /* lookup table used in tree construction */
35 double tci[8][3] = { {-1, -1, -1}, {-1, 1, -1}, {1, -1, -1}, {1, 1, -1},
36 {-1, -1, 1}, {-1, 1, 1}, {1, -1, 1}, {1, 1, 1},
37 };
38 
42 static inline int cellInNode(int c, struct bht_node *node)
43 {
44 
45  double xmin, xmax, ymin, ymax, zmin, zmax;
46  double l;
47  double x, y, z;
48  int an;
49  l = node->len * 0.5;
50  x = cells[c].x;
51  y = cells[c].y;
52  z = cells[c].z;
53 
54  xmin = node->center[0] - l;
55  ymin = node->center[1] - l;
56  if (sdim == 3)
57  zmin = node->center[2] - l;
58  if (sdim == 2)
59  zmin = 0.0;
60  xmax = node->center[0] + l;
61  ymax = node->center[1] + l;
62  if (sdim == 3)
63  zmax = node->center[2] + l;
64  if (sdim == 2)
65  zmax = 0.0;
66 
67  if (x <= xmax && x > xmin && y <= ymax && y > ymin
68  && (sdim == 2 || (z <= zmax && z > zmin)))
69  return 1;
70 
71  return 0;
72 }
73 
74 
78 static inline int whichChildNode(int p, struct bht_node *node)
79 {
80 
81  int i, j, k;
82  int num = 0;
83  double x_min, x_max, y_min, y_max, z_min, z_max;
84  double x, y, z;
85  double cx, cy, cz;
86  double shift;
87  int idx;
88 
89  shift = node->len * 0.25;
90  x = cells[p].x;
91  y = cells[p].y;
92  z = cells[p].z;
93  for (i = 0; i < tnc; i++) {
94  cx = node->center[0] + tci[i][0] * shift;
95  cy = node->center[1] + tci[i][1] * shift;
96  if (sdim == 3)
97  cz = node->center[2] + tci[i][2] * shift;
98  if (sdim == 2)
99  cz = 0.0;
100  x_min = cx - shift;
101  x_max = cx + shift;
102  y_min = cy - shift;
103  y_max = cy + shift;
104  if (sdim == 3)
105  z_min = cz - shift;
106  else
107  z_min = 0.0;
108  if (sdim == 3)
109  z_max = cz + shift;
110  else
111  z_max = 0.0;
112  if (sdim == 3) {
113  if (x > x_min && x <= x_max && y > y_min && y <= y_max && z > z_min
114  && z <= z_max)
115  break;
116  } else if (x > x_min && x <= x_max && y > y_min && y <= y_max)
117  break;
118  num++;
119  }
120 
121  if (node->child[num] == NULL) {
122 #pragma omp critical
123  {
124  idx = ni;
125  ni++;
126  }
127  node->child[num] = &tree[idx];
128  node->child[num]->len = node->len * 0.5;
129  node->child[num]->center[0] = cx;
130  node->child[num]->center[1] = cy;
131  node->child[num]->center[2] = cz;
132  node->child[num]->father = node;
133  node->child[num]->isempty = 1;
134  node->child[num]->leaf = 1;
135  for (j = 0; j < tnc; j++)
136  node->child[num]->child[j] = NULL;
137  }
138 
139  return num;
140 }
141 
145 void treeBuildInitial(int depth, struct bht_node **iroots,
146  struct bht_node *node)
147 {
148  int i, j, k;
149  if (depth != 0) {
150  double shift;
151  shift = node->len * 0.25;
152  for (j = 0; j < tnc; j++) {
153  node->child[j] = &tree[nnodes_init];
154  node->child[j]->len = node->len * 0.5;
155  node->child[j]->center[0] = node->center[0] + tci[j][0] * shift;
156  node->child[j]->center[1] = node->center[1] + tci[j][1] * shift;
157  if (sdim == 3)
158  node->child[j]->center[2] = node->center[2] + tci[j][2] * shift;
159  if (sdim == 2)
160  node->child[j]->center[2] = 0.0;
161  node->child[j]->father = node;
162  node->child[j]->isempty = 1;
163  node->child[j]->leaf = 1;
164  node->leaf = 0;
165  for (k = 0; k < tnc; k++)
166  node->child[j]->child[k] = NULL;
167  nnodes_init++;
168  treeBuildInitial(depth - 1, iroots, node->child[j]);
169  }
170  } else {
171  node->leaf = 1;
172  node->partnum = -1; /* indicates an empty node */
173  iroots[ni] = node;
174  ni++;
175  }
176 }
177 
181 void treeBuild()
182 {
183 
184  int i, j;
185  int pow2, idepth, pparts, k, l;
186  int maxTreeSize;
187  struct bht_node *node;
188  struct bht_node **iroots;
189  int nnodes_per_thread = 0;
190 
191  if (lnc == 0)
192  return;
193 #ifdef DEBUG
194  if (MPIrank == 0 && !(step % statOutStep)) {
195  printf(" Tree build...");
196  fflush(stdout);
197  }
198 #endif
199  if (sdim == 3)
200  tnc = 8;
201  if (sdim == 2)
202  tnc = 4;
203 
204  nnodes = 0;
205  nnodes_init = 0;
206 
207  pow2 = 0;
208  k = 1;
209  while (1) {
210  if (k >= OMPthreads || k >= lnc)
211  break;
212  k = k * 2;
213  pow2 += 1;
214  }
215 
216  idepth = 0;
217  l = 1;
218  while (1) {
219  if (l >= k)
220  break;
221  l = l * tnc;
222  idepth += 1;
223  }
224 
225  if (!
226  (iroots =
227  (struct bht_node **) malloc(l * sizeof(struct bht_node *))))
228  stopRun(106, "iroots", __FILE__, __LINE__);
229 
230  pparts = l / k;
231 
232  /* allocate memory for the leafs table */
233  if (!
234  (leafs =
235  (struct bht_node **) malloc(lnc * sizeof(struct bht_node *))))
236  stopRun(106, "leafs", __FILE__, __LINE__);
237 
238  maxTreeSize = lnc * 16;
239 
240  if (!
241  (tree =
242  (struct bht_node *) malloc(maxTreeSize * sizeof(struct bht_node))))
243  stopRun(106, "tree", __FILE__, __LINE__);
244 
245  bht_root = &tree[0];
246 
247  nnodes_init = 1;
248 
249 
250  boxmin[0] = nx;
251  boxmin[1] = nx;
252  boxmin[2] = nx;
253  boxmax[0] = -nx;
254  boxmax[1] = -nx;
255  boxmax[2] = -nx;
256 
257  for (i = 0; i < lnc; i++) {
258  boxmin[0] = (cells[i].x < boxmin[0] ? cells[i].x : boxmin[0]);
259  boxmax[0] = (cells[i].x > boxmax[0] ? cells[i].x : boxmax[0]);
260  boxmin[1] = (cells[i].y < boxmin[1] ? cells[i].y : boxmin[1]);
261  boxmax[1] = (cells[i].y > boxmax[1] ? cells[i].y : boxmax[1]);
262  boxmin[2] = (cells[i].z < boxmin[2] ? cells[i].z : boxmin[2]);
263  boxmax[2] = (cells[i].z > boxmax[2] ? cells[i].z : boxmax[2]);
264  }
265 
266  boxsize = boxmax[0] - boxmin[0];
267  boxsize =
268  (boxmax[1] - boxmin[1] > boxsize ? boxmax[1] - boxmin[1] : boxsize);
269  boxsize =
270  (boxmax[2] - boxmin[2] > boxsize ? boxmax[2] - boxmin[2] : boxsize);
271  boxsize += 0.1;
272 
273  /* set up the root node */
274  bht_root->len = boxsize;
275  bht_root->center[0] = boxmin[0] + (boxmax[0] - boxmin[0]) / 2.0;
276  bht_root->center[1] = boxmin[1] + (boxmax[1] - boxmin[1]) / 2.0;
277  if (sdim == 3)
278  bht_root->center[2] = boxmin[2] + (boxmax[2] - boxmin[2]) / 2.0;
279  if (sdim == 2)
280  bht_root->center[2] = 0.0;
281 
282  bht_root->leaf = 1;
283  bht_root->isempty = 1;
284  bht_root->father = NULL;
285 
286  for (i = 0; i < tnc; i++)
287  bht_root->child[i] = NULL;
288 
289  ni = 0;
290 
291  treeBuildInitial(idepth, iroots, bht_root);
292 
293  ni = nnodes_init;
294 
295  /* stop if there are no particles in this process */
296  if (lnc == 0)
297  return;
298 
299 #pragma omp parallel num_threads(k) private(node,i,j) firstprivate(pparts) shared(ni)
300  {
301  int OMPid;
302 
303  OMPid = omp_get_thread_num();
304 
305  node = bht_root;
306 
307  /* loop over cells - each cell is inserted into its own octree location */
308  for (i = 0; i < lnc; i++) {
309  int mycell = 0;
310  int c;
311  int s;
312  struct bht_node *no;
313 
314  no = node;
315 
316  for (j = 0; j < pparts; j++) {
317  mycell = cellInNode(i, iroots[OMPid * pparts + j]);
318  if (mycell == 1) {
319  no = iroots[OMPid * pparts + j];
320  break;
321  }
322  }
323 
324  if (mycell == 0)
325  continue;
326 
327  /* this loop iterates until given cell p is inserted into its location */
328  while (1) {
329  /* if this node is not a leaf we continue searching */
330  if (no->leaf == 0) {
331  c = whichChildNode(i, no);
332  no = no->child[c];
333  /* if this node is a leaf and it is not empty we need to move cells */
334  } else if (no->leaf == 1 && no->isempty == 0) {
335  for (s = 0; s < tnc; s++)
336  no->child[s] = NULL;
337  no->leaf = 0;
338  c = whichChildNode(no->partnum, no);
339  no->child[c]->partnum = no->partnum;
340  no->child[c]->isempty = 0;
341  leafs[no->partnum] = no->child[c];
342  c = whichChildNode(i, no);
343  no = no->child[c];
344  /* if this node is a leaf and it is empty we insert cell into the node */
345  } else {
346  no->partnum = i;
347  no->isempty = 0;
348  leafs[i] = no;
349  break;
350  }
351  }
352 
353  }
354 
355  } /* end of OpenMP parallel region */
356 
357  free(iroots);
358 #ifdef DEBUG
359  if (MPIrank == 0 && !(step % statOutStep)) {
360  printf("done\n");
361  fflush(stdout);
362  }
363 #endif
364 }
365 
369 void treeFree()
370 {
371  free(tree);
372  free(leafs);
373 }
int statOutStep
Definition: global.h:214
int MPIrank
Definition: global.h:134
#define lnc
Definition: global.h:102
float s
Definition: global.h:186
HYPRE_SStructVector x
Definition: tempf.c:41
void treeBuild()
Definition: tree.c:181
double tci[8][3]
Definition: tree.c:35
unsigned char leaf
Definition: global.h:284
int nnodes
Definition: global.h:292
void treeFree()
Definition: tree.c:369
contains various Timothy inline functions
contains the most important global variables, arrays and defines
int ni
Definition: global.h:297
unsigned char isempty
Definition: global.h:285
int sdim
Definition: global.h:160
int OMPthreads
Definition: global.h:137
int tnc
Definition: global.h:295
int step
Definition: global.h:173
struct bht_node * child[8]
Definition: global.h:283
double x
Definition: global.h:69
double center[3]
Definition: global.h:280
int partnum
Definition: global.h:281
struct bht_node * bht_root
Definition: global.h:289
struct bht_node ** leafs
Definition: global.h:290
struct cellData * cells
Definition: global.h:82
struct bht_node * father
Definition: global.h:282
double z
Definition: global.h:71
double boxmax[3]
Definition: global.h:244
double y
Definition: global.h:70
struct bht_node * tree
Definition: global.h:288
double boxsize
Definition: global.h:245
int nx
Definition: global.h:162
void stopRun(int ierr, char *name, char *file, int line)
Definition: utils.c:72
double len
Definition: global.h:279
void treeBuildInitial(int depth, struct bht_node **iroots, struct bht_node *node)
Definition: tree.c:145
int nnodes_init
Definition: global.h:293
double boxmin[3]
Definition: global.h:244