#include <stdlib.h>

#include "optitrust.h"

#include <stdio.h>

typedef struct {
  double x;
  double y;
  double z;
} vect;

typedef struct {
  float posX;
  float posY;
  float posZ;
  double speedX;
  double speedY;
  double speedZ;
} particle;

vect vect_add(vect v1, vect v2) {
  return {(v1.x + v2.x), (v1.y + v2.y), (v1.z + v2.z)};
}

vect vect_mul(double d, vect v) { return {(d * v.x), (d * v.y), (d * v.z)}; }

const int CHUNK_SIZE = 128;

typedef struct chunk {
  chunk *next;
  int size;
  particle items[CHUNK_SIZE];
} chunk;

typedef struct {
  chunk *front;
  chunk *back;
} bag;

chunk *atomic_read(chunk **p) {
  chunk *value;
  value = (*p);
  return value;
}

chunk *chunk_alloc() { return (chunk *)malloc(sizeof(chunk)); }

void chunk_free(chunk *c) { free(c); }

void bag_init(bag *b, int id_bag, int id_cell) {
  chunk *c = chunk_alloc();
  (c->size) = 0;
  (c->next) = NULL;
  (b->front) = c;
  (b->back) = c;
}

void bag_append(bag *b, bag *other, int id_bag, int id_cell) {
  if ((other->front)) {
    ((b->back)->next) = (other->front);
    (b->back) = (other->back);
    bag_init(other, id_bag, id_cell);
  }
}

void bag_nullify(bag *b) {
  (b->front) = NULL;
  (b->back) = NULL;
}

int bag_size(bag *b) {
  chunk *c = (b->front);
  int size = 0;
  while (c) {
    size += (c->size);
    c = (c->next);
  }
  return size;
}

void bag_add_front_chunk(bag *b) {
  chunk *c = chunk_alloc();
  (c->size) = 0;
  (c->next) = (b->front);
  (b->front) = c;
}

void bag_swap(bag *b1, bag *b2) {
  bag temp = (*b1);
  (*b1) = (*b2);
  (*b2) = temp;
}

typedef struct bag_iter {
  chunk *iter_chunk;
  int size;
  int index;
} bag_iter;

void bag_iter_load_chunk(bag_iter *it, chunk *c) {
  (it->iter_chunk) = c;
  (it->size) = (c->size);
  (it->index) = 0;
}

void bag_iter_init(bag_iter *it, bag *b) {
  bag_iter_load_chunk(it, (b->front));
}

bag_iter *bag_iter_begin(bag *b) {
  bag_iter *it = new bag_iter;
  bag_iter_init(it, b);
  return it;
}

particle *bag_iter_get(bag_iter *it) {
  return (&((it->iter_chunk)->items)[(it->index)]);
}

chunk *bag_iter_get_chunk(bag_iter *it) { return (it->iter_chunk); }

chunk *chunk_next(chunk *c, bool destructive) {
  chunk *cnext = (c->next);
  if (destructive) {
    chunk_free(c);
  }
  return cnext;
}

particle *bag_iter_next(bag_iter *it, bool destructive) {
  (it->index)++;
  if (((it->index) == (it->size))) {
    chunk *c = (it->iter_chunk);
    chunk *cnext = chunk_next(c, destructive);
    if ((cnext == NULL)) {
      return NULL;
    }
    bag_iter_load_chunk(it, cnext);
  }
  return bag_iter_get(it);
}

void bag_ho_iter_basic(bag *b, void body(particle *)) {
  bag_iter *it = bag_iter_begin(b);
  for (particle *p = bag_iter_get(it); (p != NULL);
       p = bag_iter_next(it, true)) {
    body(p);
  }
  free(it);
}

void bag_ho_iter_chunk(bag *b, void body(particle *)) {
  for (chunk *c = (b->front); (c != NULL); c = chunk_next(c, true)) {
    int nb = (c->size);
    for (int i = 0; (i < nb); i++) {
      particle *p = (&(c->items)[i]);
      body(p);
    }
  }
}

void bag_init_initial(bag *b) { bag_init(b, (-1), (-1)); }

unsigned int FREELIST_SIZE;

int **free_index;

chunk ***free_chunks;

chunk **all_free_chunks;

int number_of_spare_chunks_per_parity;

int **spare_chunks_ids;

int num_threads;

int *cumulative_free_indexes;

int bag_last_spare_chunk_to_be_used;

chunk *manual_chunk_alloc(int thread_id) {
  if ((free_index[thread_id][0] > 0)) {
    return free_chunks[thread_id][--free_index[thread_id][0]];
  } else {
    return (chunk *)malloc(sizeof(chunk));
  }
}

void manual_chunk_free(chunk *c, int thread_id) {
  if ((free_index[thread_id][0] < FREELIST_SIZE)) {
    free_chunks[thread_id][free_index[thread_id][0]++] = c;
  } else {
    free(c);
  }
}

const int THREAD_INITIAL = (-1);

const int THREAD_ZERO = 0;

chunk *manual_obtain_chunk_initial() {
  if ((free_index[THREAD_ZERO][0] < 1)) {
    fprintf(stderr,
            "Not enough chunks in all_free_chunks. Check its allocation.\n");
    exit(1);
  }
  return all_free_chunks[--free_index[THREAD_ZERO][0]];
}

chunk *manual_obtain_chunk(int id_bag, int id_cell, int thread_id) {
  if ((thread_id == THREAD_INITIAL)) {
    return manual_obtain_chunk_initial();
  }
  int id_chunk = spare_chunks_ids[id_bag][id_cell];
  int k;
  for (k = (num_threads - 1); (k >= 0); k--)
    if ((cumulative_free_indexes[k] <= id_chunk))
      break;
  int offset = (id_chunk - cumulative_free_indexes[k]);
  if (((offset < 0) || (offset >= free_index[k][0]))) {
    printf("Not enough free chunks in thread %d !\n", k);
    printf("Maybe did you forgot to call compute_cumulative_free_list_sizes "
           "and/or update_free_list_sizes ?\n");
    exit(1);
  }
  chunk *c = free_chunks[k][((free_index[k][0] - 1) - offset)];
  return c;
}

void compute_cumulative_free_list_sizes() {
  int k;
  cumulative_free_indexes[0] = 0;
  for (k = 1; (k < num_threads); k++)
    cumulative_free_indexes[k] =
        (cumulative_free_indexes[(k - 1)] + free_index[(k - 1)][0]);
  int nb_free_chunks = (cumulative_free_indexes[(num_threads - 1)] +
                        free_index[(num_threads - 1)][0]);
  if ((nb_free_chunks < number_of_spare_chunks_per_parity)) {
    int nb_chunks_to_allocate =
        (number_of_spare_chunks_per_parity - nb_free_chunks);
    printf("Not enough free chunks in the free lists ! We must malloc %d "
           "chunks.\n",
           nb_chunks_to_allocate);
    int nb_allocated_chunks = 0;
    while ((nb_allocated_chunks < nb_chunks_to_allocate)) {
      free_chunks[THREAD_ZERO][free_index[THREAD_ZERO][0]++] =
          (chunk *)malloc(sizeof(chunk));
      nb_allocated_chunks++;
    }
    compute_cumulative_free_list_sizes();
  }
  for (k = (num_threads - 1); (k >= 0); k--)
    if ((cumulative_free_indexes[k] <= number_of_spare_chunks_per_parity)) {
      bag_last_spare_chunk_to_be_used = k;
      return;
    }
}

void update_free_list_sizes() {
  for (int i = 0; (i < bag_last_spare_chunk_to_be_used); i++)
    free_index[i][0] = 0;
  free_index[bag_last_spare_chunk_to_be_used][0] -=
      (number_of_spare_chunks_per_parity -
       cumulative_free_indexes[bag_last_spare_chunk_to_be_used]);
}

bag *CHOOSE(int nb, bag *b1, bag *b2) { return b1; }

const double areaX = 10.;

const double areaY = 10.;

const double areaZ = 10.;

const double stepDuration = 0.2;

const double particleCharge = 10.;

const double particleMass = 5.;

const int gridX = 64;

const int gridY = 64;

const int gridZ = 64;

int block = 2;

int halfBlock = (block / 2);

const int nbCells = ((gridX * gridY) * gridZ);

const double cellX = (areaX / gridX);

const double cellY = (areaY / gridY);

const double cellZ = (areaZ / gridZ);

const int nbSteps = 100;

int int_of_double(double a) { return ((int)a - (a < 0.)); }

int wrap(int gridSize, int a) {
  return (((a % gridSize) + gridSize) % gridSize);
}

const int nbCorners = 8;

vect *fields = (vect *)malloc((nbCells * sizeof(vect)));

int cellOfCoord(int i, int j, int k) {
  return MINDEX3(gridX, gridY, gridZ, i, j, k);
}

int idCellOfPos(vect pos) {
  int iX = int_of_double((pos.x / cellX));
  int iY = int_of_double((pos.y / cellY));
  int iZ = int_of_double((pos.z / cellZ));
  return cellOfCoord(iX, iY, iZ);
}

double relativePosX(double x) {
  int iX = int_of_double((x / cellX));
  return ((x - (iX * cellX)) / cellX);
}

double relativePosY(double y) {
  int iY = int_of_double((y / cellY));
  return ((y - (iY * cellY)) / cellY);
}

double relativePosZ(double z) {
  int iZ = int_of_double((z / cellZ));
  return ((z - (iZ * cellZ)) / cellZ);
}

typedef struct {
  int iX;
  int iY;
  int iZ;
} coord;

coord coordOfCell(int idCell) {
  const int iZ = (idCell % gridZ);
  const int iXY = (idCell / gridZ);
  const int iY = (iXY % gridY);
  const int iX = (iXY / gridY);
  return {iX, iY, iZ};
}

typedef struct { int v[nbCorners]; } int_nbCorners;

typedef struct { double v[nbCorners]; } double_nbCorners;

typedef struct { vect v[nbCorners]; } vect_nbCorners;

int_nbCorners indicesOfCorners(int idCell) {
  const coord coord = coordOfCell(idCell);
  const int x = coord.iX;
  const int y = coord.iY;
  const int z = coord.iZ;
  const int x2 = wrap(gridX, (x + 1));
  const int y2 = wrap(gridY, (y + 1));
  const int z2 = wrap(gridZ, (z + 1));
  return {cellOfCoord(x, y, z),   cellOfCoord(x, y, z2),
          cellOfCoord(x, y2, z),  cellOfCoord(x, y2, z2),
          cellOfCoord(x2, y, z),  cellOfCoord(x2, y, z2),
          cellOfCoord(x2, y2, z), cellOfCoord(x2, y2, z2)};
}

vect_nbCorners getFieldAtCorners(int idCell, vect *field) {
  const int_nbCorners indices = indicesOfCorners(idCell);
  vect_nbCorners res;
  for (int k = 0; (k < nbCorners); k++) {
    res.v[k] = field[indices.v[k]];
  }
  return res;
}

void accumulateChargeAtCorners(double *nextCharge, int idCell,
                               double_nbCorners charges) {
  const int_nbCorners indices = indicesOfCorners(idCell);
  for (int k = 0; (k < nbCorners); k++) {
    nextCharge[indices.v[k]] += charges.v[k];
  }
}

double_nbCorners cornerInterpolationCoeff(vect pos) {
  double coefX[8] = {1., 1., 1., 1., 0., 0., 0., 0.};
  double signX[8] = {(-1.), (-1.), (-1.), (-1.), 1., 1., 1., 1.};
  double coefY[8] = {1., 1., 0., 0., 1., 1., 0., 0.};
  double signY[8] = {(-1.), (-1.), 1., 1., (-1.), (-1.), 1., 1.};
  double coefZ[8] = {1., 0., 1., 0., 1., 0., 1., 0.};
  double signZ[8] = {(-1.), 1., (-1.), 1., (-1.), 1., (-1.), 1.};
  int iX = int_of_double((pos.x / cellX));
  double rX = ((pos.x - (iX * cellX)) / cellX);
  int iY = int_of_double((pos.y / cellY));
  double rY = ((pos.y - (iY * cellY)) / cellY);
  int iZ = int_of_double((pos.z / cellZ));
  double rZ = ((pos.z - (iZ * cellZ)) / cellZ);
  double_nbCorners r;
  for (int k = 0; (k < nbCorners); k++) {
    r.v[k] = (((coefX[k] + (signX[k] * rX)) * (coefY[k] + (signY[k] * rY))) *
              (coefZ[k] + (signZ[k] * rZ)));
  }
  return r;
}

vect matrix_vect_mul(const double_nbCorners coeffs,
                     const vect_nbCorners matrix) {
  vect res = {0., 0., 0.};
  res.x =
      (res.x +
       ((((((((coeffs.v[0] * matrix.v[0].x) + (coeffs.v[1] * matrix.v[1].x)) +
             (coeffs.v[2] * matrix.v[2].x)) +
            (coeffs.v[3] * matrix.v[3].x)) +
           (coeffs.v[4] * matrix.v[4].x)) +
          (coeffs.v[5] * matrix.v[5].x)) +
         (coeffs.v[6] * matrix.v[6].x)) +
        (coeffs.v[7] * matrix.v[7].x)));
  res.y =
      (res.y +
       ((((((((coeffs.v[0] * matrix.v[0].y) + (coeffs.v[1] * matrix.v[1].y)) +
             (coeffs.v[2] * matrix.v[2].y)) +
            (coeffs.v[3] * matrix.v[3].y)) +
           (coeffs.v[4] * matrix.v[4].y)) +
          (coeffs.v[5] * matrix.v[5].y)) +
         (coeffs.v[6] * matrix.v[6].y)) +
        (coeffs.v[7] * matrix.v[7].y)));
  res.z =
      (res.z +
       ((((((((coeffs.v[0] * matrix.v[0].z) + (coeffs.v[1] * matrix.v[1].z)) +
             (coeffs.v[2] * matrix.v[2].z)) +
            (coeffs.v[3] * matrix.v[3].z)) +
           (coeffs.v[4] * matrix.v[4].z)) +
          (coeffs.v[5] * matrix.v[5].z)) +
         (coeffs.v[6] * matrix.v[6].z)) +
        (coeffs.v[7] * matrix.v[7].z)));
  return res;
}

double_nbCorners vect8_mul(const double a, const double_nbCorners data) {
  double_nbCorners res;
  for (int k = 0; (k < nbCorners); k++) {
    res.v[k] = (a * data.v[k]);
  }
  return res;
}

void init(bag *bagsCur, bag *bagsNext, vect *field);

void updateFieldUsingNextCharge(double *nextCharge, vect *field) {}

double coefX0[8] = {1., 1., 1., 1., 0., 0., 0., 0.};

double signX0[8] = {(-1.), (-1.), (-1.), (-1.), 1., 1., 1., 1.};

double coefY0[8] = {1., 1., 0., 0., 1., 1., 0., 0.};

double signY0[8] = {(-1.), (-1.), 1., 1., (-1.), (-1.), 1., 1.};

double coefZ0[8] = {1., 0., 1., 0., 1., 0., 1., 0.};

double signZ0[8] = {(-1.), 1., (-1.), 1., (-1.), 1., (-1.), 1.};

void bag_push_concurrent(bag *b, particle p) {
  chunk *c;
  int index;
  while (true) {
    c = (b->front);
    index = (c->size)++;
    if ((index < CHUNK_SIZE)) {
      (c->items)[index].posX = ((p.posX / (1. / cellX)) * (1. / cellX));
      (c->items)[index].posY = ((p.posY / (1. / cellY)) * (1. / cellY));
      (c->items)[index].posZ = ((p.posZ / (1. / cellZ)) * (1. / cellZ));
      (c->items)[index].speedX =
          ((p.speedX / (stepDuration / cellX)) * (stepDuration / cellX));
      (c->items)[index].speedY =
          ((p.speedY / (stepDuration / cellY)) * (stepDuration / cellY));
      (c->items)[index].speedZ =
          ((p.speedZ / (stepDuration / cellZ)) * (stepDuration / cellZ));
      if ((index == (CHUNK_SIZE - 1))) {
        bag_add_front_chunk(b);
      }
      return;
    } else {
      (c->size) = CHUNK_SIZE;
      while ((atomic_read((&(b->front))) == c)) {
      }
    }
  }
}

void bag_push_serial(bag *b, particle p) {
  chunk *c = (b->front);
  int index = (c->size)++;
  (c->items)[index].posX = ((p.posX / (1. / cellX)) * (1. / cellX));
  (c->items)[index].posY = ((p.posY / (1. / cellY)) * (1. / cellY));
  (c->items)[index].posZ = ((p.posZ / (1. / cellZ)) * (1. / cellZ));
  (c->items)[index].speedX =
      ((p.speedX / (stepDuration / cellX)) * (stepDuration / cellX));
  (c->items)[index].speedY =
      ((p.speedY / (stepDuration / cellY)) * (stepDuration / cellY));
  (c->items)[index].speedZ =
      ((p.speedZ / (stepDuration / cellZ)) * (stepDuration / cellZ));
  if ((index == (CHUNK_SIZE - 1))) {
    bag_add_front_chunk(b);
  }
}

void bag_push(bag *b, particle p) { return bag_push_serial(b, p); }

void bag_push_initial(bag *b, particle p) { bag_push_serial(b, p); }

int mybij(int nbCells, int nbCorners, int idCell, int idCorner) {
  coord coord = coordOfCell(idCell);
  int iX = coord.iX;
  int iY = coord.iY;
  int iZ = coord.iZ;
  int res[] = {cellOfCoord(iX, iY, iZ),
               cellOfCoord(iX, iY, wrap(gridZ, (iZ - 1))),
               cellOfCoord(iX, wrap(gridY, (iY - 1)), iZ),
               cellOfCoord(iX, wrap(gridY, (iY - 1)), wrap(gridZ, (iZ - 1))),
               cellOfCoord(wrap(gridX, (iX - 1)), iY, iZ),
               cellOfCoord(wrap(gridX, (iX - 1)), iY, wrap(gridZ, (iZ - 1))),
               cellOfCoord(wrap(gridX, (iX - 1)), wrap(gridY, (iY - 1)), iZ),
               cellOfCoord(wrap(gridX, (iX - 1)), wrap(gridY, (iY - 1)),
                           wrap(gridZ, (iZ - 1)))};
  return MINDEX2(nbCells, nbCorners, res[idCorner], idCorner);
}

int omp_get_thread_num();

int nbThreads = 8;

int main() {
  bag *bagsCur = (bag *)malloc((nbCells * sizeof(bag)));
  bag *bagsNext = (bag *)malloc((nbCells * sizeof(bag)));
  double *nextCharge = (double *)MMALLOC1(nbCells, sizeof(double));
  vect *field = (vect *)malloc((nbCells * sizeof(vect)));
  init(bagsCur, bagsNext, field);
  double *nextChargeCorners =
      (double *)MMALLOC2(nbCells, nbCorners, sizeof(double));
  double *nextChargeThreadCorners =
      (double *)MMALLOC3(nbThreads, nbCells, nbCorners, sizeof(double));
  for (int step = 0; (step < nbSteps); step++) {
    updateFieldUsingNextCharge(nextCharge, field);
    for (int idCell = 0; (idCell < nbCells); idCell++) {
      for (int idCorner = 0; (idCorner < nbCorners); idCorner++) {
        for (int k = 0; (k < nbThreads); k++) {
          nextChargeThreadCorners[MINDEX3(nbThreads, nbCells, nbCorners, k,
                                          idCell, idCorner)] = 0.;
        }
      }
    }
  core:
    for (int ciX = 0; (ciX < block); ciX++) {
      for (int ciY = 0; (ciY < block); ciY++) {
        for (int ciZ = 0; (ciZ < block); ciZ++) {
#pragma omp parallel for shared(bX, bY, bZ)
          for (int biX = (ciX * block); (biX < gridX); biX += (block * block)) {
            for (int biY = (ciY * block); (biY < gridY);
                 biY += (block * block)) {
              for (int biZ = (ciZ * block); (biZ < gridZ);
                   biZ += (block * block)) {
                for (int iX = biX; (iX < (biX + block)); iX++) {
                  for (int iY = biY; (iY < (biY + block)); iY++) {
                    for (int iZ = biZ; (iZ < (biZ + block)); iZ++) {
                      int idCell = ((((iX * gridY) + iY) * gridZ) + iZ);
                      vect_nbCorners field_at_corners =
                          getFieldAtCorners(idCell, field);
                      bag *b = (&bagsCur[idCell]);
                      for (chunk *c = (b->front); (c != NULL);
                           c = chunk_next(c, true)) {
                        int nb = (c->size);
                        int idCell2_step[nb];
                        for (int i = 0; (i < nb); i++) {
                          int iX0 = int_of_double((c->items)[i].posX);
                          double rX0 = ((c->items)[i].posX - iX0);
                          int iY0 = int_of_double((c->items)[i].posY);
                          double rY0 = ((c->items)[i].posY - iY0);
                          int iZ0 = int_of_double((c->items)[i].posZ);
                          double rZ0 = ((c->items)[i].posZ - iZ0);
                          double_nbCorners coeffs;
                          for (int k = 0; (k < nbCorners); k++) {
                            coeffs.v[k] = (((coefX0[k] + (signX0[k] * rX0)) *
                                            (coefY0[k] + (signY0[k] * rY0))) *
                                           (coefZ0[k] + (signZ0[k] * rZ0)));
                          }
                          double fieldAtPosX = 0.;
                          double fieldAtPosY = 0.;
                          double fieldAtPosZ = 0.;
                          fieldAtPosX =
                              (fieldAtPosX +
                               ((((((((coeffs.v[0] * field_at_corners.v[0].x) +
                                      (coeffs.v[1] * field_at_corners.v[1].x)) +
                                     (coeffs.v[2] * field_at_corners.v[2].x)) +
                                    (coeffs.v[3] * field_at_corners.v[3].x)) +
                                   (coeffs.v[4] * field_at_corners.v[4].x)) +
                                  (coeffs.v[5] * field_at_corners.v[5].x)) +
                                 (coeffs.v[6] * field_at_corners.v[6].x)) +
                                (coeffs.v[7] * field_at_corners.v[7].x)));
                          fieldAtPosY =
                              (fieldAtPosY +
                               ((((((((coeffs.v[0] * field_at_corners.v[0].y) +
                                      (coeffs.v[1] * field_at_corners.v[1].y)) +
                                     (coeffs.v[2] * field_at_corners.v[2].y)) +
                                    (coeffs.v[3] * field_at_corners.v[3].y)) +
                                   (coeffs.v[4] * field_at_corners.v[4].y)) +
                                  (coeffs.v[5] * field_at_corners.v[5].y)) +
                                 (coeffs.v[6] * field_at_corners.v[6].y)) +
                                (coeffs.v[7] * field_at_corners.v[7].y)));
                          fieldAtPosZ =
                              (fieldAtPosZ +
                               ((((((((coeffs.v[0] * field_at_corners.v[0].z) +
                                      (coeffs.v[1] * field_at_corners.v[1].z)) +
                                     (coeffs.v[2] * field_at_corners.v[2].z)) +
                                    (coeffs.v[3] * field_at_corners.v[3].z)) +
                                   (coeffs.v[4] * field_at_corners.v[4].z)) +
                                  (coeffs.v[5] * field_at_corners.v[5].z)) +
                                 (coeffs.v[6] * field_at_corners.v[6].z)) +
                                (coeffs.v[7] * field_at_corners.v[7].z)));
                          (c->items)[i].speedX =
                              ((c->items)[i].speedX + fieldAtPosX);
                          (c->items)[i].speedY =
                              ((c->items)[i].speedY + fieldAtPosY);
                          (c->items)[i].speedZ =
                              ((c->items)[i].speedZ + fieldAtPosZ);
                        }
                        for (int i = 0; (i < nb); i++) {
                          auto pX = (((c->items)[i].posX + iX) +
                                     (c->items)[i].speedX);
                          auto pY = (((c->items)[i].posY + iY) +
                                     (c->items)[i].speedY);
                          auto pZ = (((c->items)[i].posZ + iZ) +
                                     (c->items)[i].speedZ);
                          int iX2 = int_of_double(pX);
                          int iY2 = int_of_double(pY);
                          int iZ2 = int_of_double(pZ);
                          (c->items)[i].posX = (float)(pX - iX2);
                          (c->items)[i].posY = (float)(pY - iY2);
                          (c->items)[i].posZ = (float)(pZ - iZ2);
                          int &idCell2 = idCell2_step[i];
                          idCell2 = cellOfCoord(iX2, iY2, iZ2);
                        }
                        for (int i = 0; (i < nb); i++) {
                          int &idCell2 = (idCell2_step[i]);
                          particle p2;
                          p2.posX = (c->items)[i].posX;
                          p2.posY = (c->items)[i].posY;
                          p2.posZ = (c->items)[i].posZ;
                          p2.speedX = (c->items)[i].speedX;
                          p2.speedY = (c->items)[i].speedY;
                          p2.speedZ = (c->items)[i].speedZ;
                          const coord co = (coordOfCell(idCell2));
                          const bool isDistFromBlockLessThanHalfABlock =
                              (((co.iX - biX >= -halfBlock &&
                                 co.iX - biX < block + halfBlock) &&
                                (co.iY - biY >= -halfBlock &&
                                 co.iY - biY < block + halfBlock)) &&
                               (co.iZ - biZ >= -halfBlock &&
                                co.iZ - biZ < block + halfBlock));
                          bag *b2 = (&bagsNext[idCell2]);
                          if (isDistFromBlockLessThanHalfABlock) {
                            bag_push_serial(b2, p2);
                          } else {
                            bag_push_concurrent(b2, p2);
                          }
                          double rX1 = (c->items)[i].posX;
                          double rY1 = (c->items)[i].posY;
                          double rZ1 = (c->items)[i].posZ;
                          double_nbCorners coeffs2;
                          double_nbCorners deltaChargeOnCorners;
                          int idThread = omp_get_thread_num();
                          ;
                        charge:
                          for (int k = 0; (k < nbCorners); k++) {
                            nextChargeThreadCorners[MINDEX3(nbThreads, nbCells,
                                                            nbCorners, idThread,
                                                            idCell2, k)] +=
                                (particleCharge *
                                 (((coefX0[k] + (signX0[k] * rX1)) *
                                   (coefY0[k] + (signY0[k] * rY1))) *
                                  (coefZ0[k] + (signZ0[k] * rZ1))));
                          }
                        }
                      }
                      bag_init_initial(b);
                    }
                  }
                }
              }
            }
          }
        }
      }
    }
#pragma omp parallel for shared(idCell)
    for (int idCell = 0; (idCell < nbCells); idCell++) {
      for (int idCorner = 0; (idCorner < nbCorners); idCorner++) {
        int sum = 0;
        for (int k = 0; (k < nbThreads); k++) {
          sum += nextChargeThreadCorners[MINDEX3(nbThreads, nbCells, nbCorners,
                                                 k, idCell, idCorner)];
        }
        nextChargeCorners[MINDEX2(nbCells, nbCorners, idCell, idCorner)] = sum;
      }
    }
#pragma omp parallel for shared(idCell)
    for (int idCell = 0; (idCell < nbCells); idCell++) {
      int sum = 0;
      for (int k = 0; (k < nbCorners); k++) {
        sum += nextChargeCorners[mybij(nbCells, nbCorners, idCell, k)];
      }
      nextCharge[MINDEX1(nbCells, idCell)] = sum;
    }
    for (int idCell = 0; (idCell < nbCells); idCell++) {
      bag_swap((&bagsCur[idCell]), (&bagsNext[idCell]));
    }
  }
  MFREE(nextChargeCorners);
  MFREE(nextChargeThreadCorners);
}