cell.c 16.4 KB
Newer Older
1
/*******************************************************************************
2
 * This file is part of SWIFT.
3
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
4
 *
5
6
7
8
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
9
 *
10
11
12
13
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
14
 *
15
16
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
17
 *
18
19
20
21
22
23
24
25
26
 ******************************************************************************/

/* Config parameters. */
#include "../config.h"

/* Some standard headers. */
#include <float.h>
#include <limits.h>
#include <math.h>
27
28
29
30
#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
31

32
33
/* MPI headers. */
#ifdef WITH_MPI
34
#include <mpi.h>
35
36
#endif

37
38
/* Switch off timers. */
#ifdef TIMER
39
#undef TIMER
40
41
#endif

42
43
44
/* This object's header. */
#include "cell.h"

45
/* Local headers. */
46
#include "atomic.h"
47
#include "error.h"
48
#include "hydro.h"
49
50
#include "space.h"
#include "timers.h"
51

52
53
54
/* Global variables. */
int cell_next_tag = 0;

55
56
57
58
59
60
/**
 * @brief Get the size of the cell subtree.
 *
 * @param c The #cell.
 */

61
int cell_getsize(struct cell *c) {
62

Pedro Gonnet's avatar
Pedro Gonnet committed
63
64
  /* Number of cells in this subtree. */
  int count = 1;
65

66
67
  /* Sum up the progeny if split. */
  if (c->split)
Pedro Gonnet's avatar
Pedro Gonnet committed
68
    for (int k = 0; k < 8; k++)
69
70
71
72
73
74
75
      if (c->progeny[k] != NULL) count += cell_getsize(c->progeny[k]);

  /* Return the final count. */
  return count;
}

/**
76
77
78
79
80
81
82
83
84
 * @brief Unpack the data of a given cell and its sub-cells.
 *
 * @param pc An array of packed #pcell.
 * @param c The #cell in which to unpack the #pcell.
 * @param s The #space in which the cells are created.
 *
 * @return The number of cells created.
 */

85
86
87
88
int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {

  /* Unpack the current pcell. */
  c->h_max = pc->h_max;
89
90
  c->ti_end_min = pc->ti_end_min;
  c->ti_end_max = pc->ti_end_max;
91
  c->count = pc->count;
92
  c->gcount = pc->gcount;
93
  c->tag = pc->tag;
Matthieu Schaller's avatar
Matthieu Schaller committed
94

95
96
  /* Number of new cells created. */
  int count = 1;
97
98

  /* Fill the progeny recursively, depth-first. */
Pedro Gonnet's avatar
Pedro Gonnet committed
99
  for (int k = 0; k < 8; k++)
100
    if (pc->progeny[k] >= 0) {
Pedro Gonnet's avatar
Pedro Gonnet committed
101
      struct cell *temp = space_getcell(s);
102
      temp->count = 0;
103
      temp->gcount = 0;
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
      temp->loc[0] = c->loc[0];
      temp->loc[1] = c->loc[1];
      temp->loc[2] = c->loc[2];
      temp->h[0] = c->h[0] / 2;
      temp->h[1] = c->h[1] / 2;
      temp->h[2] = c->h[2] / 2;
      temp->dmin = c->dmin / 2;
      if (k & 4) temp->loc[0] += temp->h[0];
      if (k & 2) temp->loc[1] += temp->h[1];
      if (k & 1) temp->loc[2] += temp->h[2];
      temp->depth = c->depth + 1;
      temp->split = 0;
      temp->dx_max = 0.0;
      temp->nodeID = c->nodeID;
      temp->parent = c;
      c->progeny[k] = temp;
      c->split = 1;
      count += cell_unpack(&pc[pc->progeny[k]], temp, s);
122
123
    }

124
125
126
  /* Return the total number of unpacked cells. */
  return count;
}
127

128
/**
129
 * @brief Link the cells recursively to the given #part array.
130
131
132
133
134
135
136
 *
 * @param c The #cell.
 * @param parts The #part array.
 *
 * @return The number of particles linked.
 */

137
int cell_link_parts(struct cell *c, struct part *parts) {
138

139
140
141
  c->parts = parts;

  /* Fill the progeny recursively, depth-first. */
Pedro Gonnet's avatar
Pedro Gonnet committed
142
143
144
145
  if (c->split) {
    int offset = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL)
146
        offset += cell_link_parts(c->progeny[k], &parts[offset]);
Pedro Gonnet's avatar
Pedro Gonnet committed
147
148
    }
  }
149

150
  /* Return the total number of linked particles. */
151
152
  return c->count;
}
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
/**
 * @brief Link the cells recursively to the given #gpart array.
 *
 * @param c The #cell.
 * @param gparts The #gpart array.
 *
 * @return The number of particles linked.
 */

int cell_link_gparts(struct cell *c, struct gpart *gparts) {

  c->gparts = gparts;

  /* Fill the progeny recursively, depth-first. */
  if (c->split) {
    int offset = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL)
        offset += cell_link_gparts(c->progeny[k], &gparts[offset]);
    }
  }

  /* Return the total number of linked particles. */
  return c->gcount;
}

180
181
182
183
184
185
186
187
188
189
/**
 * @brief Pack the data of the given cell and all it's sub-cells.
 *
 * @param c The #cell.
 * @param pc Pointer to an array of packed cells in which the
 *      cells will be packed.
 *
 * @return The number of packed cells.
 */

190
191
192
193
int cell_pack(struct cell *c, struct pcell *pc) {

  /* Start by packing the data of the current cell. */
  pc->h_max = c->h_max;
194
195
  pc->ti_end_min = c->ti_end_min;
  pc->ti_end_max = c->ti_end_max;
196
  pc->count = c->count;
197
  pc->gcount = c->gcount;
198
199
200
  c->tag = pc->tag = atomic_inc(&cell_next_tag) % cell_max_tag;

  /* Fill in the progeny, depth-first recursion. */
Pedro Gonnet's avatar
Pedro Gonnet committed
201
202
  int count = 1;
  for (int k = 0; k < 8; k++)
203
204
205
206
207
208
209
210
211
    if (c->progeny[k] != NULL) {
      pc->progeny[k] = count;
      count += cell_pack(c->progeny[k], &pc[count]);
    } else
      pc->progeny[k] = -1;

  /* Return the number of packed cells used. */
  return count;
}
212

213
214
215
216
217
218
/**
 * @brief Lock a cell and hold its parents.
 *
 * @param c The #cell.
 */

219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
int cell_locktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to lock this cell. */
  if (c->hold || lock_trylock(&c->lock) != 0) {
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
  if (c->hold) {

    /* Unlock this cell. */
    if (lock_unlock(&c->lock) != 0) error("Failed to unlock cell.");

    /* Admit defeat. */
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Climb up the tree and lock/hold/unlock. */
Pedro Gonnet's avatar
Pedro Gonnet committed
241
  struct cell *finger;
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
  for (finger = c->parent; finger != NULL; finger = finger->parent) {

    /* Lock this cell. */
    if (lock_trylock(&finger->lock) != 0) break;

    /* Increment the hold. */
    atomic_inc(&finger->hold);

    /* Unlock the cell. */
    if (lock_unlock(&finger->lock) != 0) error("Failed to unlock cell.");
  }

  /* If we reached the top of the tree, we're done. */
  if (finger == NULL) {
    TIMER_TOC(timer_locktree);
    return 0;
  }

  /* Otherwise, we hit a snag. */
  else {

    /* Undo the holds up to finger. */
Pedro Gonnet's avatar
Pedro Gonnet committed
264
265
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
      __sync_fetch_and_sub(&finger2->hold, 1);

    /* Unlock this cell. */
    if (lock_unlock(&c->lock) != 0) error("Failed to unlock cell.");

    /* Admit defeat. */
    TIMER_TOC(timer_locktree);
    return 1;
  }
}

int cell_glocktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to lock this cell. */
  if (c->ghold || lock_trylock(&c->glock) != 0) {
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
  if (c->ghold) {

    /* Unlock this cell. */
    if (lock_unlock(&c->glock) != 0) error("Failed to unlock cell.");

    /* Admit defeat. */
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Climb up the tree and lock/hold/unlock. */
Pedro Gonnet's avatar
Pedro Gonnet committed
299
  struct cell *finger;
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
  for (finger = c->parent; finger != NULL; finger = finger->parent) {

    /* Lock this cell. */
    if (lock_trylock(&finger->glock) != 0) break;

    /* Increment the hold. */
    __sync_fetch_and_add(&finger->ghold, 1);

    /* Unlock the cell. */
    if (lock_unlock(&finger->glock) != 0) error("Failed to unlock cell.");
  }

  /* If we reached the top of the tree, we're done. */
  if (finger == NULL) {
    TIMER_TOC(timer_locktree);
    return 0;
  }

  /* Otherwise, we hit a snag. */
  else {

    /* Undo the holds up to finger. */
Pedro Gonnet's avatar
Pedro Gonnet committed
322
323
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
324
325
326
327
328
329
330
331
332
333
      __sync_fetch_and_sub(&finger2->ghold, 1);

    /* Unlock this cell. */
    if (lock_unlock(&c->glock) != 0) error("Failed to unlock cell.");

    /* Admit defeat. */
    TIMER_TOC(timer_locktree);
    return 1;
  }
}
334

335
/**
336
 * @brief Unlock a cell's parents.
337
338
339
 *
 * @param c The #cell.
 */
340
341
342
343
344
345
346
347
348

void cell_unlocktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to unlock this cell. */
  if (lock_unlock(&c->lock) != 0) error("Failed to unlock cell.");

  /* Climb up the tree and unhold the parents. */
Pedro Gonnet's avatar
Pedro Gonnet committed
349
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
350
351
352
353
354
355
356
357
358
359
360
361
362
    __sync_fetch_and_sub(&finger->hold, 1);

  TIMER_TOC(timer_locktree);
}

void cell_gunlocktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to unlock this cell. */
  if (lock_unlock(&c->glock) != 0) error("Failed to unlock cell.");

  /* Climb up the tree and unhold the parents. */
Pedro Gonnet's avatar
Pedro Gonnet committed
363
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
364
365
366
367
368
    __sync_fetch_and_sub(&finger->ghold, 1);

  TIMER_TOC(timer_locktree);
}

369
370
371
372
373
/**
 * @brief Sort the parts into eight bins along the given pivots.
 *
 * @param c The #cell array to be sorted.
 */
374
375
376

void cell_split(struct cell *c) {

Pedro Gonnet's avatar
Pedro Gonnet committed
377
378
  int i, j;
  const int count = c->count, gcount = c->gcount;
379
380
381
  struct part *parts = c->parts;
  struct xpart *xparts = c->xparts;
  struct gpart *gparts = c->gparts;
382
383
384
385
  int left[8], right[8];
  double pivot[3];

  /* Init the pivots. */
Pedro Gonnet's avatar
Pedro Gonnet committed
386
  for (int k = 0; k < 3; k++) pivot[k] = c->loc[k] + c->h[k] / 2;
387
388
389
390
391
392
393
394

  /* Split along the x-axis. */
  i = 0;
  j = count - 1;
  while (i <= j) {
    while (i <= count - 1 && parts[i].x[0] <= pivot[0]) i += 1;
    while (j >= 0 && parts[j].x[0] > pivot[0]) j -= 1;
    if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
395
      struct part temp = parts[i];
396
397
      parts[i] = parts[j];
      parts[j] = temp;
Pedro Gonnet's avatar
Pedro Gonnet committed
398
      struct xpart xtemp = xparts[i];
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
      xparts[i] = xparts[j];
      xparts[j] = xtemp;
    }
  }
  /* for ( k = 0 ; k <= j ; k++ )
      if ( parts[k].x[0] > pivot[0] )
          error( "cell_split: sorting failed." );
  for ( k = i ; k < count ; k++ )
      if ( parts[k].x[0] < pivot[0] )
          error( "cell_split: sorting failed." ); */
  left[1] = i;
  right[1] = count - 1;
  left[0] = 0;
  right[0] = j;

  /* Split along the y axis, twice. */
Pedro Gonnet's avatar
Pedro Gonnet committed
415
  for (int k = 1; k >= 0; k--) {
416
417
418
419
420
421
    i = left[k];
    j = right[k];
    while (i <= j) {
      while (i <= right[k] && parts[i].x[1] <= pivot[1]) i += 1;
      while (j >= left[k] && parts[j].x[1] > pivot[1]) j -= 1;
      if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
422
        struct part temp = parts[i];
423
424
        parts[i] = parts[j];
        parts[j] = temp;
Pedro Gonnet's avatar
Pedro Gonnet committed
425
        struct xpart xtemp = xparts[i];
426
427
428
429
430
431
432
433
        xparts[i] = xparts[j];
        xparts[j] = xtemp;
      }
    }
    /* for ( int kk = left[k] ; kk <= j ; kk++ )
        if ( parts[kk].x[1] > pivot[1] ) {
            message( "ival=[%i,%i], i=%i, j=%i." , left[k] , right[k] , i , j );
            error( "sorting failed (left)." );
434
            }
435
436
437
438
439
440
441
442
443
444
    for ( int kk = i ; kk <= right[k] ; kk++ )
        if ( parts[kk].x[1] < pivot[1] )
            error( "sorting failed (right)." ); */
    left[2 * k + 1] = i;
    right[2 * k + 1] = right[k];
    left[2 * k] = left[k];
    right[2 * k] = j;
  }

  /* Split along the z axis, four times. */
Pedro Gonnet's avatar
Pedro Gonnet committed
445
  for (int k = 3; k >= 0; k--) {
446
447
448
449
450
451
    i = left[k];
    j = right[k];
    while (i <= j) {
      while (i <= right[k] && parts[i].x[2] <= pivot[2]) i += 1;
      while (j >= left[k] && parts[j].x[2] > pivot[2]) j -= 1;
      if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
452
        struct part temp = parts[i];
453
454
        parts[i] = parts[j];
        parts[j] = temp;
Pedro Gonnet's avatar
Pedro Gonnet committed
455
        struct xpart xtemp = xparts[i];
456
457
458
459
460
461
462
463
        xparts[i] = xparts[j];
        xparts[j] = xtemp;
      }
    }
    /* for ( int kk = left[k] ; kk <= j ; kk++ )
        if ( parts[kk].x[2] > pivot[2] ) {
            message( "ival=[%i,%i], i=%i, j=%i." , left[k] , right[k] , i , j );
            error( "sorting failed (left)." );
464
            }
465
466
467
468
469
470
471
472
473
474
475
476
    for ( int kk = i ; kk <= right[k] ; kk++ )
        if ( parts[kk].x[2] < pivot[2] ) {
            message( "ival=[%i,%i], i=%i, j=%i." , left[k] , right[k] , i , j );
            error( "sorting failed (right)." );
            } */
    left[2 * k + 1] = i;
    right[2 * k + 1] = right[k];
    left[2 * k] = left[k];
    right[2 * k] = j;
  }

  /* Store the counts and offsets. */
Pedro Gonnet's avatar
Pedro Gonnet committed
477
  for (int k = 0; k < 8; k++) {
478
479
480
481
482
483
    c->progeny[k]->count = right[k] - left[k] + 1;
    c->progeny[k]->parts = &c->parts[left[k]];
    c->progeny[k]->xparts = &c->xparts[left[k]];
  }

  /* Re-link the gparts. */
Pedro Gonnet's avatar
Pedro Gonnet committed
484
  for (int k = 0; k < count; k++)
485
486
487
488
489
490
491
492
493
494
495
496
497
    if (parts[k].gpart != NULL) parts[k].gpart->part = &parts[k];

  /* Verify that _all_ the parts have been assigned to a cell. */
  /* for ( k = 1 ; k < 8 ; k++ )
      if ( &c->progeny[k-1]->parts[ c->progeny[k-1]->count ] !=
  c->progeny[k]->parts )
          error( "Particle sorting failed (internal consistency)." );
  if ( c->progeny[0]->parts != c->parts )
      error( "Particle sorting failed (left edge)." );
  if ( &c->progeny[7]->parts[ c->progeny[7]->count ] != &c->parts[ count ] )
      error( "Particle sorting failed (right edge)." ); */

  /* Verify a few sub-cells. */
Pedro Gonnet's avatar
Pedro Gonnet committed
498
  /* for (int k = 0 ; k < c->progeny[0]->count ; k++ )
499
500
501
502
      if ( c->progeny[0]->parts[k].x[0] > pivot[0] ||
           c->progeny[0]->parts[k].x[1] > pivot[1] ||
           c->progeny[0]->parts[k].x[2] > pivot[2] )
          error( "Sorting failed (progeny=0)." );
Pedro Gonnet's avatar
Pedro Gonnet committed
503
  for (int k = 0 ; k < c->progeny[1]->count ; k++ )
504
505
506
507
      if ( c->progeny[1]->parts[k].x[0] > pivot[0] ||
           c->progeny[1]->parts[k].x[1] > pivot[1] ||
           c->progeny[1]->parts[k].x[2] <= pivot[2] )
          error( "Sorting failed (progeny=1)." );
Pedro Gonnet's avatar
Pedro Gonnet committed
508
  for (int k = 0 ; k < c->progeny[2]->count ; k++ )
509
510
511
512
513
514
515
516
517
518
519
520
521
522
      if ( c->progeny[2]->parts[k].x[0] > pivot[0] ||
           c->progeny[2]->parts[k].x[1] <= pivot[1] ||
           c->progeny[2]->parts[k].x[2] > pivot[2] )
          error( "Sorting failed (progeny=2)." ); */

  /* Now do the same song and dance for the gparts. */

  /* Split along the x-axis. */
  i = 0;
  j = gcount - 1;
  while (i <= j) {
    while (i <= gcount - 1 && gparts[i].x[0] <= pivot[0]) i += 1;
    while (j >= 0 && gparts[j].x[0] > pivot[0]) j -= 1;
    if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
523
      struct gpart gtemp = gparts[i];
524
525
      gparts[i] = gparts[j];
      gparts[j] = gtemp;
526
    }
527
528
529
530
531
532
533
  }
  left[1] = i;
  right[1] = gcount - 1;
  left[0] = 0;
  right[0] = j;

  /* Split along the y axis, twice. */
Pedro Gonnet's avatar
Pedro Gonnet committed
534
  for (int k = 1; k >= 0; k--) {
535
536
537
538
539
540
    i = left[k];
    j = right[k];
    while (i <= j) {
      while (i <= right[k] && gparts[i].x[1] <= pivot[1]) i += 1;
      while (j >= left[k] && gparts[j].x[1] > pivot[1]) j -= 1;
      if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
541
        struct gpart gtemp = gparts[i];
542
543
544
545
546
547
548
549
550
551
552
        gparts[i] = gparts[j];
        gparts[j] = gtemp;
      }
    }
    left[2 * k + 1] = i;
    right[2 * k + 1] = right[k];
    left[2 * k] = left[k];
    right[2 * k] = j;
  }

  /* Split along the z axis, four times. */
Pedro Gonnet's avatar
Pedro Gonnet committed
553
  for (int k = 3; k >= 0; k--) {
554
555
556
557
558
559
    i = left[k];
    j = right[k];
    while (i <= j) {
      while (i <= right[k] && gparts[i].x[2] <= pivot[2]) i += 1;
      while (j >= left[k] && gparts[j].x[2] > pivot[2]) j -= 1;
      if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
560
        struct gpart gtemp = gparts[i];
561
562
563
564
565
566
567
568
569
570
571
        gparts[i] = gparts[j];
        gparts[j] = gtemp;
      }
    }
    left[2 * k + 1] = i;
    right[2 * k + 1] = right[k];
    left[2 * k] = left[k];
    right[2 * k] = j;
  }

  /* Store the counts and offsets. */
Pedro Gonnet's avatar
Pedro Gonnet committed
572
  for (int k = 0; k < 8; k++) {
573
574
575
576
577
    c->progeny[k]->gcount = right[k] - left[k] + 1;
    c->progeny[k]->gparts = &c->gparts[left[k]];
  }

  /* Re-link the parts. */
Pedro Gonnet's avatar
Pedro Gonnet committed
578
  for (int k = 0; k < gcount; k++)
579
580
    if (gparts[k].id > 0) gparts[k].part->gpart = &gparts[k];
}
581
582
583
584
585
586
587
588
589
590
591

/**
 * @brief Initialises all particles to a valid state even if the ICs were stupid
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
void cell_init_parts(struct cell *c, void *data) {

  struct part *p = c->parts;
  struct xpart *xp = c->xparts;
592
  const int count = c->count;
593

594
  for (int i = 0; i < count; ++i) {
595
596
    p[i].ti_begin = 0;
    p[i].ti_end = 0;
597
598
599
    xp[i].v_full[0] = p[i].v[0];
    xp[i].v_full[1] = p[i].v[1];
    xp[i].v_full[2] = p[i].v[2];
600
    hydro_first_init_part(&p[i], &xp[i]);
601
602
    hydro_init_part(&p[i]);
    hydro_reset_acceleration(&p[i]);
603
  }
604
605
  c->ti_end_min = 0;
  c->ti_end_max = 0;
606
607
}

608
/**
609
610
 * @brief Converts hydro quantities to a valid state after the initial density
 *calculation
611
612
613
614
615
616
617
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
void cell_convert_hydro(struct cell *c, void *data) {

  struct part *p = c->parts;
618
619

  for (int i = 0; i < c->count; ++i) {
620
621
622
623
    hydro_convert_quantities(&p[i]);
  }
}

Matthieu Schaller's avatar
Matthieu Schaller committed
624
625
626
627
628
629
/**
 * @brief Cleans the links in a given cell.
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
630
void cell_clean_links(struct cell *c, void *data) {
Matthieu Schaller's avatar
Matthieu Schaller committed
631
632
  c->density = NULL;
  c->nr_density = 0;
633

Matthieu Schaller's avatar
Matthieu Schaller committed
634
635
636
  c->force = NULL;
  c->nr_force = 0;
}