ParticleFactory.cpph

Go to the documentation of this file.

#include <random>
 
#include "peano4/utils/Loop.h"
 
 
namespace {
  std::random_device               rd;
  std::mt19937                     e2(rd());
  std::uniform_real_distribution<> randomDist(0, 1);
} // namespace
 
 
template <class T>
void toolbox::particles::init(T& particle, const tarch::la::Vector<Dimensions, double>& x, double searchRadius) {
#if PeanoDebug > 0
  particle.setDebugX(x);
  particle.setDebugH(0.0);
#endif
 
  particle.setX(x);
  particle.setSearchRadius(searchRadius);
 
  particle.setParallelState(T::ParallelState::Virtual);
}
 
 
template <class T>
std::vector<T*> toolbox::particles::createEquallySpacedParticles(
  double                                       spacing,
  const tarch::la::Vector<Dimensions, double>& voxelX,
  const tarch::la::Vector<Dimensions, double>& voxelH,
  bool                                         addNoise
) {
  std::vector<T*> result;
 
  double continuousParticlesPerAxis = voxelH(0) / spacing;
  int    particlesPerAxis           = static_cast<int>(std::floor(continuousParticlesPerAxis));
  continuousParticlesPerAxis -= particlesPerAxis;
  assertion4(continuousParticlesPerAxis >= 0.0, continuousParticlesPerAxis, particlesPerAxis, voxelH(0), spacing);
  assertion4(continuousParticlesPerAxis <= 1.0, continuousParticlesPerAxis, particlesPerAxis, voxelH(0), spacing);
 
  const double areaOfGridSpannedByCellLocalCartesianGridAlongOneAxis = particlesPerAxis * spacing;
  assertion(tarch::la::smallerEquals(areaOfGridSpannedByCellLocalCartesianGridAlongOneAxis, voxelH(0)));
 
  if (particlesPerAxis > 0) {
    const double particleSpacing = voxelH(0) / particlesPerAxis;
 
    const tarch::la::Vector<Dimensions, double>
      offset = voxelX - 0.5 * areaOfGridSpannedByCellLocalCartesianGridAlongOneAxis;
 
    dfor(k, particlesPerAxis) {
      tarch::la::Vector<Dimensions, double> x = offset + particleSpacing * tarch::la::convertScalar<double>(k);
      T*                                    newParticle = new T();
 
      if (addNoise) {
        for (int d = 0; d < Dimensions; d++) {
          x(d) += (randomDist(e2) - 0.5) * particleSpacing;
        }
      }
 
      init(*newParticle, x, 0.9 * tarch::la::min(voxelH));
      result.push_back(newParticle);
    }
  }
 
  return result;
}
 
 
template <class T>
std::vector<T*> toolbox::particles::createParticlesAlignedWithGlobalCartesianMesh(
  double                                       spacing,
  const tarch::la::Vector<Dimensions, double>& voxelX,
  const tarch::la::Vector<Dimensions, double>& voxelH,
  const tarch::la::Vector<Dimensions, double>& domainOffset,
  const tarch::la::Vector<Dimensions, double>& domainSize,
  bool                                         addNoise
) {
  static tarch::logging::Log _log("toolbox::particles");
  logTraceInWith4Arguments(
    "createParticlesAlignedWithGlobalCartesianMesh(...)", spacing, voxelX, domainOffset, addNoise
  );
  std::vector<T*> result;
 
  tarch::la::Vector<Dimensions, double> biasedOffset = domainOffset + 0.5 * tarch::la::remainder(domainOffset, spacing);
 
  tarch::la::Vector<Dimensions, int>    particlesPerAxis;
  tarch::la::Vector<Dimensions, double> firstParticleX;
  for (int d = 0; d < Dimensions; d++) {
    int leftParticleNumberAlongAxis = static_cast<int>(
      std::round((voxelX(d) - 0.5 * voxelH(d) - biasedOffset(d)) / spacing)
    );
    int rightParticleNumberAlongAxis = static_cast<int>(
      std::round((voxelX(d) + 0.5 * voxelH(d) - biasedOffset(d)) / spacing)
    );
 
    double leftX  = biasedOffset(d) + spacing * leftParticleNumberAlongAxis;
    double rightX = biasedOffset(d) + spacing * rightParticleNumberAlongAxis;
 
    // One of them has to be biased
    if (tarch::la::smaller(leftX, voxelX(d) - 0.5 * voxelH(d)))
      leftParticleNumberAlongAxis++;
    if (tarch::la::equals(rightX, voxelX(d) + 0.5 * voxelH(d)) and tarch::la::smaller(rightX, domainOffset(d) + domainSize(d)) or tarch::la::greater(rightX, voxelX(d) + 0.5 * voxelH(d)))
      rightParticleNumberAlongAxis--;
 
    particlesPerAxis(d) = std::max(0, rightParticleNumberAlongAxis - leftParticleNumberAlongAxis + 1);
    firstParticleX(d)   = biasedOffset(d) + spacing * leftParticleNumberAlongAxis;
  }
 
  const double RelativeAccuracyForGeometricChecks = 0.01 * voxelH(0);
 
  dfor(k, particlesPerAxis) {
    tarch::la::Vector<Dimensions, double> x = firstParticleX + spacing * tarch::la::convertScalar<double>(k);
 
    T* newParticle = new T();
 
    if (addNoise) {
      for (int d = 0; d < Dimensions; d++) {
        x(d) += (randomDist(e2) - 0.5) * spacing;
      }
    }
 
    init(*newParticle, x, 0.9 * tarch::la::min(voxelH));
    result.push_back(newParticle);
  }
 
  logTraceOutWith1Argument("createParticlesAlignedWithGlobalCartesianMesh(...)", result.size());
  return result;
}