1 #include "objects/finders/UniformGrid.h"
2 #include "system/Profiler.h"
3
4 NAMESPACE_SPH_BEGIN
5
UniformGridFinder(const Float relativeCellCnt)6 UniformGridFinder::UniformGridFinder(const Float relativeCellCnt)
7 : relativeCellCnt(relativeCellCnt) {
8 Indices::init();
9 }
10
11 UniformGridFinder::~UniformGridFinder() = default;
12
buildImpl(IScheduler & UNUSED (scheduler),ArrayView<const Vector> points)13 void UniformGridFinder::buildImpl(IScheduler& UNUSED(scheduler), ArrayView<const Vector> points) {
14 PROFILE_SCOPE("VoxelFinder::buildImpl");
15 // number of voxels, free parameter
16 const Size lutSize = Size(relativeCellCnt * root<3>(points.size())) + 1;
17 if (lut.empty() || lutSize != lut.getDimensionSize()) {
18 // build lookup map if not yet build or we have a significantly different number of points
19 lut = LookupMap(lutSize);
20 }
21 if (SPH_LIKELY(!points.empty())) {
22 lut.update(points);
23 }
24 }
25
26 template <bool FindAll>
find(const Vector & pos,const Size index,const Float radius,Array<NeighborRecord> & neighbors) const27 Size UniformGridFinder::find(const Vector& pos,
28 const Size index,
29 const Float radius,
30 Array<NeighborRecord>& neighbors) const {
31 const Vector refPosition = lut.clamp(pos);
32 Indices lower = lut.map(refPosition);
33 Indices upper = lower;
34 Box voxel = lut.voxel(lower);
35 const Vector size = lut.getVoxelSize();
36 Vector diffUpper = voxel.upper() - pos;
37 Vector diffLower = pos - voxel.lower();
38
39 SPH_ASSERT(lut.getDimensionSize() > 0);
40 const int upperLimit = lut.getDimensionSize() - 1;
41 while (upper[X] < upperLimit && diffUpper[X] < radius) {
42 diffUpper[X] += size[X];
43 upper[X]++;
44 }
45 while (lower[X] > 0 && diffLower[X] < radius) {
46 diffLower[X] += size[X];
47 lower[X]--;
48 }
49 while (upper[Y] < upperLimit && diffUpper[Y] < radius) {
50 diffUpper[Y] += size[Y];
51 upper[Y]++;
52 }
53 while (lower[Y] > 0 && diffLower[Y] < radius) {
54 diffLower[Y] += size[Y];
55 lower[Y]--;
56 }
57 while (upper[Z] < upperLimit && diffUpper[Z] < radius) {
58 diffUpper[Z] += size[Z];
59 upper[Z]++;
60 }
61 while (lower[Z] > 0 && diffLower[Z] < radius) {
62 diffLower[Z] += size[Z];
63 lower[Z]--;
64 }
65
66 for (int x = lower[X]; x <= upper[X]; ++x) {
67 for (int y = lower[Y]; y <= upper[Y]; ++y) {
68 for (int z = lower[Z]; z <= upper[Z]; ++z) {
69 for (Size i : lut(Indices(x, y, z))) {
70 const Float distSqr = getSqrLength(values[i] - pos);
71 if (distSqr < sqr(radius) && (FindAll || rank[i] < rank[index])) {
72 neighbors.emplaceBack(NeighborRecord{ i, distSqr });
73 }
74 }
75 }
76 }
77 }
78 return neighbors.size();
79 }
80
81 template Size UniformGridFinder::find<true>(const Vector& pos,
82 const Size index,
83 const Float radius,
84 Array<NeighborRecord>& neighbors) const;
85
86 template Size UniformGridFinder::find<false>(const Vector& pos,
87 const Size index,
88 const Float radius,
89 Array<NeighborRecord>& neighbors) const;
90
91 NAMESPACE_SPH_END
92