CellBlock3D.h

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
#include <array>
#include <cmath>
#include <numeric>
#include <vector>
 
#include "autopas/containers/CellBorderAndFlagManager.h"
#include "autopas/options/IteratorBehavior.h"
#include "autopas/utils/ArrayMath.h"
#include "autopas/utils/ExceptionHandler.h"
#include "autopas/utils/ThreeDimensionalMapping.h"
#include "autopas/utils/inBox.h"
 
namespace autopas::internal {
template <class ParticleCell>
class CellBlock3D : public CellBorderAndFlagManager {
 public:
  using index_t = std::size_t;
  CellBlock3D(std::vector<ParticleCell> &vec, const std::array<double, 3> &bMin, const std::array<double, 3> &bMax,
              double interactionLength, double cellSizeFactor = 1.0)
      : _cells(&vec), _boxMin(bMin), _boxMax(bMax), _interactionLength(interactionLength) {
    rebuild(vec, bMin, bMax, interactionLength, cellSizeFactor);
 
    for (int i = 0; i < 3; ++i) {
      if (bMax[i] < bMin[i] + interactionLength) {
        AutoPasLog(ERROR, "Interaction length too large is {}, bmin {}, bmax {}", interactionLength, bMin[i], bMax[i]);
        utils::ExceptionHandler::exception("Error in CellBlock3D: interaction Length too large!");
      }
    }
  }
 
  CellBlock3D(const CellBlock3D &) = delete;
 
  CellBlock3D &operator=(const CellBlock3D) = delete;
 
  [[nodiscard]] bool cellCanContainHaloParticles(index_t index1d) const override {
    if (index1d < _firstOwnedCellIndex or index1d > _lastOwnedCellIndex) {
      return true;
    }
    const auto index3d = oneToThreeD(index1d);
    bool isHaloCell = false;
    for (size_t i = 0; i < index3d.size(); ++i) {
      if (index3d[i] < _cellsPerInteractionLength or
          index3d[i] >= _cellsPerDimensionWithHalo[i] - _cellsPerInteractionLength) {
        isHaloCell = true;
        break;
      }
    }
    return isHaloCell;
  }
 
  [[nodiscard]] bool cellCanContainOwnedParticles(index_t index1d) const override {
    return not cellCanContainHaloParticles(index1d);
  }
 
  ParticleCell &getCell(index_t index1d) const;
 
  ParticleCell &getCell(const std::array<index_t, 3> &index3d) const;
 
  void reserve(size_t numParticles);
 
  void rebuild(std::vector<ParticleCell> &vec, const std::array<double, 3> &bMin, const std::array<double, 3> &bMax,
               double interactionLength, double cellSizeFactor);
 
  ParticleCell &getContainingCell(const std::array<double, 3> &pos) const;
 
  std::pair<std::array<double, 3>, std::array<double, 3>> getCellBoundingBox(index_t index1d) const;
 
  std::pair<std::array<double, 3>, std::array<double, 3>> getCellBoundingBox(
      const std::array<index_t, 3> &index3d) const;
 
  [[nodiscard]] std::array<index_t, 3> get3DIndexOfPosition(const std::array<double, 3> &pos) const;
 
  [[nodiscard]] index_t get1DIndexOfPosition(const std::array<double, 3> &pos) const;
 
  [[nodiscard]] const std::array<index_t, 3> &getCellsPerDimensionWithHalo() const {
    return _cellsPerDimensionWithHalo;
  }
 
  [[nodiscard]] bool checkInHalo(const std::array<double, 3> &position) const;
 
  void clearHaloCells();
 
  std::vector<ParticleCell *> getNearbyHaloCells(const std::array<double, 3> &position, double allowedDistance) const {
    using namespace autopas::utils::ArrayMath::literals;
 
    const auto index3D = get3DIndexOfPosition(position);
    const auto index1D = threeToOneD(index3D);
 
    // these indices are (already) at least 0 and at most _cellsPerDimensionWithHalo[i]-1
    const auto lowIndex3D = get3DIndexOfPosition(position - allowedDistance);
    const auto highIndex3D = get3DIndexOfPosition(position + allowedDistance);
 
    std::vector<ParticleCell *> closeHaloCells;
    // This is an overestimation with the upper bound of possible number of cells in the vicinity.
    // An overestimate is cheaper than many reallocations.
    // If the size of the overallocation ever becomes a problem we can use vector::shrink_to_fit() before return.
    const auto blockLength = highIndex3D - lowIndex3D;
    const auto numInterestingCells = std::max(1ul, blockLength[0] * blockLength[1] * blockLength[2]);
    closeHaloCells.reserve(numInterestingCells);
 
    // always add the cell the particle is currently in first if it is a halo cell.
    if ((*_cells)[index1D].getPossibleParticleOwnerships() == OwnershipState::halo) {
      closeHaloCells.push_back(&getCell(index3D));
    }
 
    auto currentIndex3D = lowIndex3D;
    for (currentIndex3D[0] = lowIndex3D[0]; currentIndex3D[0] <= highIndex3D[0]; ++currentIndex3D[0]) {
      for (currentIndex3D[1] = lowIndex3D[1]; currentIndex3D[1] <= highIndex3D[1]; ++currentIndex3D[1]) {
        for (currentIndex3D[2] = lowIndex3D[2]; currentIndex3D[2] <= highIndex3D[2]; ++currentIndex3D[2]) {
          // we have already added the cell which normally would belong to the particle, so skip here:
          if (currentIndex3D == index3D) {
            continue;
          }
          // we need to return the cell is it is a halo cell.
          const auto currentIndex1D = threeToOneD(currentIndex3D);
          if ((*_cells)[currentIndex1D].getPossibleParticleOwnerships() == OwnershipState::halo) {
            closeHaloCells.push_back(&getCell(currentIndex3D));
          }
        }
      }
    }
 
    return closeHaloCells;
  }
 
  [[nodiscard]] const std::array<double, 3> &getHaloBoxMin() const { return _haloBoxMin; }
 
  [[nodiscard]] const std::array<double, 3> &getHaloBoxMax() const { return _haloBoxMax; }
 
  [[nodiscard]] unsigned long getCellsPerInteractionLength() const { return _cellsPerInteractionLength; }
 
  [[nodiscard]] const std::array<double, 3> &getCellLength() const { return _cellLength; }
 
  index_t getNumCells() const;
 
  index_t getFirstOwnedCellIndex() const;
 
  index_t getLastOwnedCellIndex() const;
 
 private:
  [[nodiscard]] std::array<index_t, 3> oneToThreeD(index_t index1D) const;
 
  [[nodiscard]] index_t threeToOneD(const std::array<index_t, 3> &index3D) const;
 
  std::array<index_t, 3> _cellsPerDimensionWithHalo{};
  index_t _firstOwnedCellIndex{};
  index_t _lastOwnedCellIndex{};
  std::vector<ParticleCell> *_cells;
 
  std::array<double, 3> _boxMin, _boxMax;
  std::array<double, 3> _haloBoxMin{}, _haloBoxMax{};
 
  double _interactionLength{};
 
  unsigned long _cellsPerInteractionLength{};
 
  std::array<double, 3> _cellLength{};
 
  std::array<double, 3> _cellLengthReciprocal{};
};
 
template <class ParticleCell>
inline typename CellBlock3D<ParticleCell>::index_t CellBlock3D<ParticleCell>::get1DIndexOfPosition(
    const std::array<double, 3> &pos) const {
  return threeToOneD(get3DIndexOfPosition(pos));
}
 
template <class ParticleCell>
inline std::array<typename CellBlock3D<ParticleCell>::index_t, 3> CellBlock3D<ParticleCell>::get3DIndexOfPosition(
    const std::array<double, 3> &pos) const {
  std::array<typename CellBlock3D<ParticleCell>::index_t, 3> cellIndex{};
 
  for (size_t dim = 0; dim < 3; dim++) {
    // 0 <= cellIndex < cellsPerDimWithHalo
    const long int value = (static_cast<long int>(std::floor((pos[dim] - _boxMin[dim]) * _cellLengthReciprocal[dim]))) +
                           _cellsPerInteractionLength;
    const index_t nonNegativeValue = static_cast<index_t>(std::max(value, 0l));
    const index_t nonLargerValue = std::min(nonNegativeValue, _cellsPerDimensionWithHalo[dim] - 1);
    cellIndex[dim] = nonLargerValue;
 
    // sanity checks and precautions for numerical instabilities around the edge of the container box,
    // which would lead to doubling of particles
    if (pos[dim] >= _boxMax[dim]) {
      // pos is located outside the box. Make sure that cellIndex is at least in the positive halo range.
      cellIndex[dim] = std::max(cellIndex[dim], _cellsPerDimensionWithHalo[dim] - _cellsPerInteractionLength);
    } else if (pos[dim] < _boxMin[dim] and cellIndex[dim] == _cellsPerInteractionLength) {
      // pos is located outside but the index resolved to be the last IN the box. Fix by decrementing the index.
      --cellIndex[dim];
    } else if (pos[dim] < _boxMax[dim] and
               cellIndex[dim] == _cellsPerDimensionWithHalo[dim] - _cellsPerInteractionLength) {
      // pos is located inside the box, but cellIndex resolves to be the first halo cell. Fix by decrementing the index.
      --cellIndex[dim];
    }
  }
 
  return cellIndex;
}
 
template <class ParticleCell>
void CellBlock3D<ParticleCell>::reserve(size_t numParticles) {
  const auto particlesPerCell = numParticles / _cells->size();
  for (auto &cell : *_cells) {
    cell.reserve(particlesPerCell);
  }
}
 
template <class ParticleCell>
inline void CellBlock3D<ParticleCell>::rebuild(std::vector<ParticleCell> &vec, const std::array<double, 3> &bMin,
                                               const std::array<double, 3> &bMax, double interactionLength,
                                               double cellSizeFactor) {
  using namespace autopas::utils::ArrayMath::literals;
 
  _cells = &vec;
  _boxMin = bMin;
  _boxMax = bMax;
  _interactionLength = interactionLength;
 
  if (cellSizeFactor >= 1.0) {
    _cellsPerInteractionLength = 1;
  } else {
    _cellsPerInteractionLength = ceil(1.0 / cellSizeFactor);
  }
  // compute cell length and number of cells
  index_t numCells = 1;
  for (int d = 0; d < 3; ++d) {
    const double boxLength = _boxMax[d] - _boxMin[d];
    // The number of cells is rounded down because the cells will be stretched to fit.
    // std::max to ensure there is at least one cell.
    const auto cellsPerDim =
        std::max(static_cast<index_t>(std::floor(boxLength / (_interactionLength * cellSizeFactor))), 1ul);
 
    _cellsPerDimensionWithHalo[d] = cellsPerDim + 2 * _cellsPerInteractionLength;
 
    _cellLength[d] = boxLength / static_cast<double>(cellsPerDim);
 
    _cellLengthReciprocal[d] = static_cast<double>(cellsPerDim) / boxLength;  // compute with least rounding possible
 
    _haloBoxMin[d] = _boxMin[d] - _cellsPerInteractionLength * _cellLength[d];
    _haloBoxMax[d] = _boxMax[d] + _cellsPerInteractionLength * _cellLength[d];
 
    numCells *= _cellsPerDimensionWithHalo[d];
  }
  AutoPasLog(TRACE, "Box Length incl Halo : {}", autopas::utils::ArrayUtils::to_string(_haloBoxMax - _haloBoxMin));
  AutoPasLog(TRACE, "Cells/Dim  incl Halo : {}", autopas::utils::ArrayUtils::to_string(_cellsPerDimensionWithHalo));
  AutoPasLog(TRACE, "Cell Length          : {}", autopas::utils::ArrayUtils::to_string(_cellLength));
  AutoPasLog(TRACE, "Interaction Length   : {}", interactionLength);
  AutoPasLog(TRACE, "Cell Size Factor     : {}", cellSizeFactor);
 
  _firstOwnedCellIndex = _cellsPerDimensionWithHalo[0] * _cellsPerDimensionWithHalo[1] * _cellsPerInteractionLength +
                         _cellsPerDimensionWithHalo[0] * _cellsPerInteractionLength + _cellsPerInteractionLength;
  _lastOwnedCellIndex = numCells - 1 - _firstOwnedCellIndex;
 
  // initialize cells
  _cells->resize(numCells);
 
  for (auto &cell : *_cells) {
    cell.setCellLength(_cellLength);
  }
 
  // determine the OwnershipStates, each cell can contain. This is later used in the CellFunctor to skip calculations
  for (int i = 0; i < numCells; i++) {
    const bool canHaveHalos = cellCanContainHaloParticles(i);
    const bool canHaveOwned = cellCanContainOwnedParticles(i);
    if (canHaveHalos and canHaveOwned) {
      (*_cells)[i].setPossibleParticleOwnerships(OwnershipState::owned | OwnershipState::halo);
    } else if (canHaveHalos) {
      (*_cells)[i].setPossibleParticleOwnerships(OwnershipState::halo);
    } else if (canHaveOwned) {
      (*_cells)[i].setPossibleParticleOwnerships(OwnershipState::owned);
    } else {
      (*_cells)[i].setPossibleParticleOwnerships(OwnershipState::dummy);
    }
  }
}
 
template <class ParticleCell>
typename CellBlock3D<ParticleCell>::index_t CellBlock3D<ParticleCell>::getNumCells() const {
  return _cells->size();
}
 
template <class ParticleCell>
typename CellBlock3D<ParticleCell>::index_t CellBlock3D<ParticleCell>::getFirstOwnedCellIndex() const {
  return _firstOwnedCellIndex;
}
 
template <class ParticleCell>
typename CellBlock3D<ParticleCell>::index_t CellBlock3D<ParticleCell>::getLastOwnedCellIndex() const {
  return _lastOwnedCellIndex;
}
 
template <class ParticleCell>
inline std::pair<std::array<double, 3>, std::array<double, 3>> CellBlock3D<ParticleCell>::getCellBoundingBox(
    const index_t index1d) const {
  return this->getCellBoundingBox(this->oneToThreeD(index1d));
}
 
template <class ParticleCell>
inline std::pair<std::array<double, 3>, std::array<double, 3>> CellBlock3D<ParticleCell>::getCellBoundingBox(
    const std::array<index_t, 3> &index3d) const {
  std::array<double, 3> boxmin{}, boxmax{};
  for (int d = 0; d < 3; d++) {
    // defaults
    boxmin[d] = index3d[d] * this->_cellLength[d] + _haloBoxMin[d];
 
    // stupid rounding errors. Snap values to the exact box values.
    if (index3d[d] == 0) {
      boxmin[d] = _haloBoxMin[d];
    } else if (index3d[d] == _cellsPerInteractionLength) {
      // Case: we are at the lower boundary of the non-halo box
      boxmin[d] = _boxMin[d];
    }
 
    boxmax[d] = boxmin[d] + this->_cellLength[d];
 
    // This must not be an else to the if block above as both cases can be true.
    // e.g. if there is only one cell
    if (index3d[d] == this->_cellsPerDimensionWithHalo[d] - _cellsPerInteractionLength - 1) {
      boxmax[d] = _boxMax[d];
    } else if (index3d[d] == this->_cellsPerDimensionWithHalo[d] - 1) {
      boxmin[d] = _haloBoxMax[d] - this->_cellLength[d];
      boxmax[d] = _haloBoxMax[d];
    }
  }
  return {boxmin, boxmax};
}
 
template <class ParticleCell>
inline ParticleCell &CellBlock3D<ParticleCell>::getContainingCell(const std::array<double, 3> &pos) const {
  auto ind = get1DIndexOfPosition(pos);
  return getCell(ind);
}
 
template <class ParticleCell>
inline ParticleCell &CellBlock3D<ParticleCell>::getCell(index_t index1d) const {
  return (*_cells)[index1d];
}
 
template <class ParticleCell>
inline ParticleCell &CellBlock3D<ParticleCell>::getCell(const std::array<index_t, 3> &index3d) const {
  return (*_cells)[threeToOneD(index3d)];
}
 
template <class ParticleCell>
inline std::array<typename CellBlock3D<ParticleCell>::index_t, 3> CellBlock3D<ParticleCell>::oneToThreeD(
    index_t index1D) const {
  return utils::ThreeDimensionalMapping::oneToThreeD(index1D, _cellsPerDimensionWithHalo);
}
 
template <class ParticleCell>
inline typename CellBlock3D<ParticleCell>::index_t CellBlock3D<ParticleCell>::threeToOneD(
    const std::array<index_t, 3> &index3D) const {
  return utils::ThreeDimensionalMapping::threeToOneD(index3D, _cellsPerDimensionWithHalo);
}
 
template <class ParticleCell>
bool CellBlock3D<ParticleCell>::checkInHalo(const std::array<double, 3> &position) const {
  return autopas::utils::inBox(position, _haloBoxMin, _haloBoxMax) &&
         autopas::utils::notInBox(position, _boxMin, _boxMax);
}
 
template <class ParticleCell>
void CellBlock3D<ParticleCell>::clearHaloCells() {
  std::vector<index_t> haloSlices(2 * _cellsPerInteractionLength);
  auto mid = haloSlices.begin() + _cellsPerInteractionLength;
  std::iota(haloSlices.begin(), mid, 0);
 
  // x: min and max of x
  std::iota(mid, haloSlices.end(), _cellsPerDimensionWithHalo[0] - _cellsPerInteractionLength);
  for (index_t i : haloSlices) {
    for (index_t j = 0; j < _cellsPerDimensionWithHalo[1]; j++) {
      for (index_t k = 0; k < _cellsPerDimensionWithHalo[2]; k++) {
        index_t index = threeToOneD({i, j, k});
        (*_cells)[index].clear();
      }
    }
  }
  // y: min and max of y
  std::iota(mid, haloSlices.end(), _cellsPerDimensionWithHalo[1] - _cellsPerInteractionLength);
  for (index_t i = 1; i < _cellsPerDimensionWithHalo[0] - 1; i++) {  // 0 and cells-1 already done in previous loop
    for (index_t j : haloSlices) {
      for (index_t k = 0; k < _cellsPerDimensionWithHalo[2]; k++) {
        index_t index = threeToOneD({i, j, k});
        (*_cells)[index].clear();
      }
    }
  }
  // z: min and max of z
  std::iota(mid, haloSlices.end(), _cellsPerDimensionWithHalo[2] - _cellsPerInteractionLength);
  for (index_t i = 1; i < _cellsPerDimensionWithHalo[0] - 1; i++) {    // 0 and cells-1 already done in previous loop
    for (index_t j = 1; j < _cellsPerDimensionWithHalo[1] - 1; j++) {  // 0 and cells-1 already done in previous loop
      for (index_t k : haloSlices) {
        index_t index = threeToOneD({i, j, k});
        (*_cells)[index].clear();
      }
    }
  }
}
}  // namespace autopas::internal