Peano 4
Loading...
Searching...
No Matches
ParticleFactory.cpph
Go to the documentation of this file.
1#include <random>
2
3#include "peano4/utils/Loop.h"
4
5
6namespace {
7 std::random_device rd;
8 std::mt19937 e2(rd());
9 std::uniform_real_distribution<> randomDist(0, 1);
10} // namespace
11
12
13template <class T>
14void toolbox::particles::init(T& particle, const tarch::la::Vector<Dimensions, double>& x, double searchRadius) {
15#if PeanoDebug > 0
16 particle.setDebugX(x);
17 particle.setDebugH(0.0);
18#endif
19
20 particle.setX(x);
21 particle.setSearchRadius(searchRadius);
22
23 particle.setParallelState(T::ParallelState::Virtual);
24}
25
26
27template <class T>
29 double spacing,
32 bool addNoise
33) {
34 std::vector<T*> result;
35
36 double continuousParticlesPerAxis = voxelH(0) / spacing;
37 int particlesPerAxis = static_cast<int>(std::floor(continuousParticlesPerAxis));
38 continuousParticlesPerAxis -= particlesPerAxis;
39 assertion4(continuousParticlesPerAxis >= 0.0, continuousParticlesPerAxis, particlesPerAxis, voxelH(0), spacing);
40 assertion4(continuousParticlesPerAxis <= 1.0, continuousParticlesPerAxis, particlesPerAxis, voxelH(0), spacing);
41
42 const double areaOfGridSpannedByCellLocalCartesianGridAlongOneAxis = particlesPerAxis * spacing;
43 assertion(tarch::la::smallerEquals(areaOfGridSpannedByCellLocalCartesianGridAlongOneAxis, voxelH(0)));
44
45 if (particlesPerAxis > 0) {
46 const double particleSpacing = voxelH(0) / particlesPerAxis;
47
49 offset = voxelX - 0.5 * areaOfGridSpannedByCellLocalCartesianGridAlongOneAxis;
50
51 dfor(k, particlesPerAxis) {
52 tarch::la::Vector<Dimensions, double> x = offset + particleSpacing * tarch::la::convertScalar<double>(k);
53 T* newParticle = new T();
54
55 if (addNoise) {
56 for (int d = 0; d < Dimensions; d++) {
57 x(d) += (randomDist(e2) - 0.5) * particleSpacing;
58 }
59 }
60
61 init(*newParticle, x, 0.9 * tarch::la::min(voxelH));
62 result.push_back(newParticle);
63 }
64 }
65
66 return result;
67}
68
69
70template <class T>
72 double spacing,
75 const tarch::la::Vector<Dimensions, double>& domainOffset,
77 bool addNoise
78) {
79 static tarch::logging::Log _log("toolbox::particles");
81 "createParticlesAlignedWithGlobalCartesianMesh(...)", spacing, voxelX, domainOffset, addNoise
82 );
83 std::vector<T*> result;
84
85 tarch::la::Vector<Dimensions, double> biasedOffset = domainOffset + 0.5 * tarch::la::remainder(domainOffset, spacing);
86
89 for (int d = 0; d < Dimensions; d++) {
90 int leftParticleNumberAlongAxis = static_cast<int>(
91 std::round((voxelX(d) - 0.5 * voxelH(d) - biasedOffset(d)) / spacing)
92 );
93 int rightParticleNumberAlongAxis = static_cast<int>(
94 std::round((voxelX(d) + 0.5 * voxelH(d) - biasedOffset(d)) / spacing)
95 );
96
97 double leftX = biasedOffset(d) + spacing * leftParticleNumberAlongAxis;
98 double rightX = biasedOffset(d) + spacing * rightParticleNumberAlongAxis;
99
100 // One of them has to be biased
101 if (tarch::la::smaller(leftX, voxelX(d) - 0.5 * voxelH(d)))
102 leftParticleNumberAlongAxis++;
103 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)))
104 rightParticleNumberAlongAxis--;
105
106 particlesPerAxis(d) = std::max(0, rightParticleNumberAlongAxis - leftParticleNumberAlongAxis + 1);
107 firstParticleX(d) = biasedOffset(d) + spacing * leftParticleNumberAlongAxis;
108 }
109
110 const double RelativeAccuracyForGeometricChecks = 0.01 * voxelH(0);
111
112 dfor(k, particlesPerAxis) {
113 tarch::la::Vector<Dimensions, double> x = firstParticleX + spacing * tarch::la::convertScalar<double>(k);
114
115 T* newParticle = new T();
116
117 if (addNoise) {
118 for (int d = 0; d < Dimensions; d++) {
119 x(d) += (randomDist(e2) - 0.5) * spacing;
120 }
121 }
122
123 init(*newParticle, x, 0.9 * tarch::la::min(voxelH));
124 result.push_back(newParticle);
125 }
126
127 logTraceOutWith1Argument("createParticlesAlignedWithGlobalCartesianMesh(...)", result.size());
128 return result;
129}
#define assertion4(expr, param0, param1, param2, param3)
#define assertion(expr)
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
#define logTraceOutWith1Argument(methodName, argument0)
Definition Log.h:380
#define logTraceInWith4Arguments(methodName, argument0, argument1, argument2, argument3)
Definition Log.h:373
#define dfor(counter, max)
d-dimensional Loop
Definition Loop.h:308
Log Device.
Definition Log.h:516
Vector< Size, Scalar > remainder(const Vector< Size, Scalar > &vector, double h)
Return the remainder of a division.
static bool smaller(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE) InlineMethod
Smaller operator for floating point values.
bool greater(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
bool smallerEquals(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
bool equals(const Matrix< Rows, Cols, Scalar > &lhs, const Matrix< Rows, Cols, Scalar > &rhs, const Scalar &tolerance=NUMERICAL_ZERO_DIFFERENCE)
Compares to matrices on equality by means of a numerical accuracy.
Scalar min(const Vector< Size, Scalar > &vector)
Returns the element with minimal value (NOT absolute value).
std::vector< T * > createParticlesAlignedWithGlobalCartesianMesh(double spacingH, 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)
Insert particles that are aligned along a Cartesian grid globally.
std::vector< T * > createEquallySpacedParticles(double spacingH, const tarch::la::Vector< Dimensions, double > &voxelX, const tarch::la::Vector< Dimensions, double > &voxelH, bool addNoise)
Create equally spaced particles for one cell.
void init(T &particle, const tarch::la::Vector< Dimensions, double > &x, double searchRadius)
Init particle.
tarch::logging::Log _log("examples::unittests")
Simple vector class.
Definition Vector.h:134