LCC18Traversal.h

Go to the documentation of this file.

 
#pragma once
 
#include "LCTraversalInterface.h"
#include "autopas/baseFunctors/CellFunctor.h"
#include "autopas/containers/cellTraversals/C18BasedTraversal.h"
#include "autopas/options/DataLayoutOption.h"
#include "autopas/utils/ArrayUtils.h"
#include "autopas/utils/ThreeDimensionalMapping.h"
#include "autopas/utils/WrapOpenMP.h"
 
namespace autopas {
 
template <class ParticleCell, class PairwiseFunctor>
class LCC18Traversal : public C18BasedTraversal<ParticleCell, PairwiseFunctor>, public LCTraversalInterface {
 public:
  explicit LCC18Traversal(const std::array<unsigned long, 3> &dims, PairwiseFunctor *pairwiseFunctor,
                          const double interactionLength, const std::array<double, 3> &cellLength,
                          DataLayoutOption dataLayout, bool useNewton3)
      : C18BasedTraversal<ParticleCell, PairwiseFunctor>(dims, pairwiseFunctor, interactionLength, cellLength,
                                                         dataLayout, useNewton3),
        _cellFunctor(pairwiseFunctor, interactionLength /*should use cutoff here, if not used to build verlet-lists*/,
                     dataLayout, useNewton3) {
    computeOffsets();
  }
 
  void traverseParticles() override;
 
  void processBaseCell(std::vector<ParticleCell> &cells, unsigned long x, unsigned long y, unsigned long z);
 
  [[nodiscard]] TraversalOption getTraversalType() const override { return TraversalOption::lc_c18; }
 
  [[nodiscard]] bool isApplicable() const override { return true; }
 
  void setSortingThreshold(size_t sortingThreshold) override { _cellFunctor.setSortingThreshold(sortingThreshold); }
 
 private:
  void computeOffsets();
 
  internal::CellFunctor<ParticleCell, PairwiseFunctor,
                        /*bidirectional*/ true>
      _cellFunctor;
 
  using offsetArray_t = std::vector<std::pair<unsigned long, std::array<double, 3>>>;
 
  std::vector<std::vector<offsetArray_t>> _cellOffsets;
 
  unsigned long getIndex(const unsigned long pos, const unsigned int dim) const;
};
 
template <class ParticleCell, class PairwiseFunctor>
inline void LCC18Traversal<ParticleCell, PairwiseFunctor>::computeOffsets() {
  _cellOffsets.resize(2 * this->_overlap[1] + 1, std::vector<offsetArray_t>(2 * this->_overlap[0] + 1));
  const std::array<long, 3> _overlap_s = utils::ArrayUtils::static_cast_copy_array<long>(this->_overlap);
 
  const auto interactionLengthSquare(this->_interactionLength * this->_interactionLength);
 
  for (long z = 0l; z <= _overlap_s[2]; ++z) {
    for (long y = -_overlap_s[1]; y <= _overlap_s[1]; ++y) {
      for (long x = -_overlap_s[0]; x <= _overlap_s[0]; ++x) {
        const long offset = utils::ThreeDimensionalMapping::threeToOneD(
            x, y, z, utils::ArrayUtils::static_cast_copy_array<long>(this->_cellsPerDimension));
 
        if (offset < 0l) {
          continue;
        }
        // add to each applicable special case
        for (long yArray = -_overlap_s[1]; yArray <= _overlap_s[1]; ++yArray) {
          if (std::abs(yArray + y) <= _overlap_s[1]) {
            for (long xArray = -_overlap_s[0]; xArray <= _overlap_s[0]; ++xArray) {
              if (std::abs(xArray + x) <= _overlap_s[0]) {
                const std::array<double, 3> pos = {
                    std::max(0l, (std::abs(x) - 1l)) * this->_cellLength[0],
                    std::max(0l, (std::abs(y) - 1l)) * this->_cellLength[1],
                    std::max(0l, (std::abs(z) - 1l)) * this->_cellLength[2],
                };
                // calculate distance between the borders of the base cell and the other cell
                const double distSquare = utils::ArrayMath::dot(pos, pos);
                // only add cell offset if cell is within cutoff radius
                if (distSquare <= interactionLengthSquare) {
                  // Calculate the sorting direction from the base cell and the other cell by use of the offset (x, y,
                  // z).
                  // Note: We have to calculate the sorting direction separately from pos, since pos is the offset
                  // between the borders of cells. For neighbouring cells this would be 0 and the sorting direction
                  // would be wrong.
                  std::array<double, 3> sortingDir = {static_cast<double>(x) * this->_cellLength[0],
                                                      static_cast<double>(y) * this->_cellLength[1],
                                                      static_cast<double>(z) * this->_cellLength[2]};
                  if (x == 0 and y == 0 and z == 0) {
                    sortingDir = {1., 1., 1.};
                  }
 
                  _cellOffsets[yArray + _overlap_s[1]][xArray + _overlap_s[0]].push_back(
                      std::make_pair(offset, utils::ArrayMath::normalize(sortingDir)));
                }
              }
            }
          }
        }
      }
    }
  }
}
 
template <class ParticleCell, class PairwiseFunctor>
unsigned long LCC18Traversal<ParticleCell, PairwiseFunctor>::getIndex(const unsigned long pos,
                                                                      const unsigned int dim) const {
  unsigned long index;
  if (pos < this->_overlap[dim]) {
    index = pos;
  } else if (pos < this->_cellsPerDimension[dim] - this->_overlap[dim]) {
    index = this->_overlap[dim];
  } else {
    index = pos - this->_cellsPerDimension[dim] + 2 * this->_overlap[dim] + 1ul;
  }
  return index;
}
 
template <class ParticleCell, class PairwiseFunctor>
void LCC18Traversal<ParticleCell, PairwiseFunctor>::processBaseCell(std::vector<ParticleCell> &cells, unsigned long x,
                                                                    unsigned long y, unsigned long z) {
  const unsigned long baseIndex = utils::ThreeDimensionalMapping::threeToOneD(x, y, z, this->_cellsPerDimension);
 
  const unsigned long xArray = getIndex(x, 0);
  const unsigned long yArray = getIndex(y, 1);
 
  ParticleCell &baseCell = cells[baseIndex];
  offsetArray_t &offsets = this->_cellOffsets[yArray][xArray];
  for (auto const &[offset, r] : offsets) {
    unsigned long otherIndex = baseIndex + offset;
    ParticleCell &otherCell = cells[otherIndex];
 
    if (baseIndex == otherIndex) {
      this->_cellFunctor.processCell(baseCell);
    } else {
      this->_cellFunctor.processCellPair(baseCell, otherCell, r);
    }
  }
}
 
template <class ParticleCell, class PairwiseFunctor>
inline void LCC18Traversal<ParticleCell, PairwiseFunctor>::traverseParticles() {
  auto &cells = *(this->_cells);
  this->template c18Traversal</*allCells*/ false>(
      [&](unsigned long x, unsigned long y, unsigned long z) { this->processBaseCell(cells, x, y, z); });
}
 
}  // namespace autopas