Skip to content

Commit 1022e29

Browse files
author
Dan Koschier
committed
Initial commit.
1 parent 56be9fe commit 1022e29

21 files changed

+2361
-2
lines changed

‎CMakeLists.txt

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,46 @@
1+
cmake_minimum_required(VERSION 3.2)
2+
3+
project(CompactNSearch)
4+
5+
# Visual studio solution directories.
6+
set_property(GLOBAL PROPERTY USE_FOLDERS on)
7+
8+
set(CMAKE_CXX_STANDARD 11)
9+
set(CMAKE_CXX_STANDARD_REQUIRED ON)
10+
11+
if (UNIX)
12+
find_package(OpenMP)
13+
if (OPENMP_FOUND)
14+
set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
15+
endif(OPENMP_FOUND)
16+
endif (UNIX)
17+
18+
SET(CMAKE_DEBUG_POSTFIX "_d")
19+
20+
set (HEADER_FILES
21+
include/Config.h
22+
include/CompactNSearch.h
23+
include/PointSet.h
24+
include/DataStructures.h)
25+
26+
include_directories("include")
27+
28+
add_library(CompactNSearch
29+
${HEADER_FILES}
30+
src/CompactNSearch.cpp
31+
)
32+
33+
install(FILES "include/CompactNSearch" ${HEADER_FILES}
34+
DESTINATION include/)
35+
36+
install(TARGETS CompactNSearch
37+
RUNTIME DESTINATION bin
38+
LIBRARY DESTINATION lib
39+
ARCHIVE DESTINATION lib)
40+
41+
option(BUILD_DEMO "Build example of how to use this library."
42+
ON)
43+
if(BUILD_DEMO)
44+
add_subdirectory(demo)
45+
endif(BUILD_DEMO)
46+

‎LICENSE

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
1-
MIT License
1+
The MIT License (MIT)
22

3-
Copyright (c) 2016 Interactive Computer Graphics
3+
Copyright (c) 2016-present, CompactNSearch contributors
44

55
Permission is hereby granted, free of charge, to any person obtaining a copy
66
of this software and associated documentation files (the "Software"), to deal
@@ -19,3 +19,4 @@ AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
1919
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
2020
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
2121
SOFTWARE.
22+

‎README.md

Lines changed: 75 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,75 @@
1+
# CompactNSearch
2+
3+
**CompactNSearch** is a C++ library for **parallel** computation of neighboring points in a **fixed radius** in a **three-dimensional point cloud**. The implemented algorithm is a variant of the *Compact Hashing* approach proposed by Ihmsen et al. [IABT11]. The neighborhood information can be efficiently updated when the points spatially move. Moreover, the library offers the possibility to reorder the points (and other array-stored per point information) according to a space-filling Z curve to improve the cache efficiency in later queries or accesses.
4+
5+
The library was used to generate all results of the SPH-based fluid simulations presented by Bender and Koschier in [BK15, BK16].
6+
7+
**Author**: Dan Koschier, **License**: MIT
8+
9+
## Libraries using CompactNSearch
10+
* [PBD] - A C++ library for physically-based simulation of rigid bodies, deformables, cloth and fluids using Position-Based Dynamics
11+
* [SPlisHSPlasH] - A C++ library for the physically-based simulation of fluids using Smoothed Particle Hydrodynamics (cf. video) (Coming soon)
12+
13+
[![Video](https://img.youtube.com/vi/POnmzzhc5E0/0.jpg)](https://www.youtube.com/watch?v=POnmzzhc5E0)
14+
15+
## Build Instructions
16+
17+
This project is based on [CMake](https://cmake.org/). Simply generate project, Makefiles, etc. using [CMake](https://cmake.org/) and compile the project with the compiler of your choice. The code was tested with the following configurations:
18+
- Windows 10 64-bit, CMake 3.7, Visual Studio 2015
19+
- Debian 8 64-bit, CMake 3.7, GCC 4.9.2.
20+
21+
## Usage
22+
A data structure to perform a neighborhood search can be created by calling the constructor given a fixed search radius ```r```.
23+
```c++
24+
CompactNSearch::NeighborhoodSearch nsearch(r);
25+
```
26+
An arbitrary number of point clouds can then be added to the data structure using the method ```add_point_set```. The library expects the point positions to be contiguously stored in an array-like structure. The method will return a unique id associated with the initialized point set.
27+
```c++
28+
std::vector<std::array<Real, 3>> positions;
29+
// ... Fill array with 3 * n real numbers representing three-dimensional point positions.
30+
unsigned int point_set_id = nsearch.add_point_set(positions.front().data(), positions.size());
31+
nsearch.find_neighbors();
32+
```
33+
In order to generate the neighborhood information simply execute the following command
34+
```c++
35+
nsearch.find_neighbors();
36+
```
37+
Finally, the neighborhood information can be accessed as follows
38+
```c++
39+
PointSet const& ps = nsearch.point_set(point_set_id);
40+
for (int i = 0; i < ps.n_points(); ++i)
41+
{
42+
for (int j = 0; j < ps.n_neighbors(i); ++j)
43+
{
44+
// Return PointID of the jth neighbor of the ith particle.
45+
PointID const& pid = ps.neighbor(i, j);
46+
// ...
47+
// Do whatever you want with the point id. The id contains two indices.
48+
// The first field pid.point_set_id represents the unique point set id returnd by add_point_set.
49+
// The second field pid.point_id stands for the index of the neighboring particle within
50+
// the containing point set.
51+
// ...
52+
}
53+
}
54+
```
55+
56+
Besides the basic functionality the library offers to compute a rule for reordering the points according to a space-filling Z curve. The reordering will improve the performance of future neighborhood queries and accesses. The rule can be computed via
57+
```c++
58+
nsearch.z_sort();
59+
```
60+
Please note that the actual reordering must invoked by the user by
61+
```c++
62+
ps.sort_field(positions.data());
63+
```
64+
Assuming that there is additional information stored per-point (e.g. velocity, color, mass etc.) the information **must** also be reorded using the same method to maintain consistency. Subsequently, the ```find_neighbors``` function has to be invoked again to update the neighborhood information.
65+
66+
Another self-explaining (benchmark) [demo](demo/main.cpp) is contained in the project.
67+
68+
## References
69+
70+
* [IABT11] M. Ihmsen, N. Akinci, M. Becker and M. Teschner, 2011. "A Parallel SPH Implementation on Multi-Core CPUs", Computer Graphics Forum 30, 1, 99-112.
71+
* [BK15] J. Bender and D. Koschier 2015. "Divergence-Free Smoothed Particle Hydrodynamics", ACM SIGGRAPH / Eurographics Symposium on Computer Animation, 1-9
72+
* [BK16] J. Bender and D. Koschier, 2016. "Divergence-Free SPH for Incompressible and Viscous Fluids", IEEE Transactions on Visualization and Computer Graphics.
73+
74+
[PBD]: <https://github.com/InteractiveComputerGraphics/PositionBasedDynamics>
75+
[SPlisHSPlasH]: <https://github.com/InteractiveComputerGraphics/SPlisHSPlasH>

‎demo/CMakeLists.txt

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
2+
3+
add_definitions(-D_USE_MATH_DEFINES)
4+
5+
include_directories("../include")
6+
7+
add_executable(Demo
8+
main.cpp
9+
)
10+
11+
target_link_libraries(Demo CompactNSearch)

‎demo/main.cpp

Lines changed: 212 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,212 @@
1+
2+
#include <CompactNSearch>
3+
4+
#include <iostream>
5+
#include <vector>
6+
#include <array>
7+
#include <cmath>
8+
#include <limits>
9+
#include <chrono>
10+
#include <random>
11+
12+
#include <omp.h>
13+
14+
using namespace CompactNSearch;
15+
16+
std::vector<std::array<Real, 3>> positions;
17+
18+
std::size_t const N = 150;
19+
Real const r_omega = 0.15;
20+
Real const r_omega2 = r_omega * r_omega;
21+
Real const radius = 2.0 * (2.0 * r_omega / static_cast<Real>(N-1));
22+
23+
std::size_t const N_enright_steps = 50;
24+
25+
Real
26+
compute_average_number_of_neighbors(NeighborhoodSearch const& nsearch)
27+
{
28+
unsigned long res = 0;
29+
auto const& d = nsearch.point_set(0);
30+
for (int i = 0; i < d.n_points(); ++i)
31+
{
32+
res += static_cast<unsigned long>(d.n_neighbors(i));
33+
}
34+
return static_cast<Real>(res) / static_cast<Real>(d.n_points());
35+
}
36+
37+
Real
38+
compute_average_distance(NeighborhoodSearch const& nsearch)
39+
{
40+
unsigned long long res = 0;
41+
auto const& d = nsearch.point_set(0);
42+
unsigned long long count = 0;
43+
for (int i = 0; i < d.n_points(); ++i)
44+
{
45+
std::size_t nn = d.n_neighbors(i);
46+
for (int j = 0; j < nn; ++j)
47+
{
48+
CompactNSearch::PointID const& k = d.neighbor(i, j);
49+
res += std::abs(i - static_cast<int>(k.point_id));
50+
count++;
51+
}
52+
}
53+
return static_cast<Real>(res) / static_cast<Real>(count);
54+
}
55+
56+
std::vector<std::vector<unsigned int>>
57+
brute_force_search()
58+
{
59+
std::vector<std::vector<unsigned int>> brute_force_neighbors(positions.size());
60+
for (int i = 0; i < positions.size(); ++i)
61+
{
62+
std::vector<unsigned int>& neighbors = brute_force_neighbors[i];
63+
for (int j = 0; j < positions.size(); ++j)
64+
{
65+
if (i == j)
66+
continue;
67+
std::array<Real, 3> const& xa = positions[i];
68+
std::array<Real, 3> const& xb = positions[j];
69+
Real l2 =
70+
(xa[0] - xb[0])*(xa[0] - xb[0]) +
71+
(xa[1] - xb[1])*(xa[1] - xb[1]) +
72+
(xa[2] - xb[2])*(xa[2] - xb[2]);
73+
if (l2 <= radius * radius)
74+
{
75+
neighbors.push_back(j);
76+
}
77+
}
78+
}
79+
return std::move(brute_force_neighbors);
80+
}
81+
82+
void
83+
compare_with_bruteforce_search(NeighborhoodSearch const& nsearch)
84+
{
85+
auto brute_force_neighbors = brute_force_search();
86+
PointSet const& d0 = nsearch.point_set(0);
87+
for (int i = 0; i < N; ++i)
88+
{
89+
auto const& bfn = brute_force_neighbors[i];
90+
if (bfn.size() != d0.n_neighbors(i))
91+
{
92+
std::cerr << "ERROR: Not the same number of neighbors." << std::endl;
93+
}
94+
for (int j = 0; j < d0.n_neighbors(i); ++j)
95+
{
96+
if (std::find(bfn.begin(), bfn.end(), d0.neighbor(i, j).point_id) == bfn.end())
97+
{
98+
std::cerr << "ERROR: Neighbor not found in brute force list." << std::endl;
99+
}
100+
}
101+
}
102+
}
103+
104+
std::array<Real, 3>
105+
enright_velocity_field(std::array<Real, 3> const& x)
106+
{
107+
Real sin_pi_x_2 = std::sin(M_PI * x[0]);
108+
Real sin_pi_y_2 = std::sin(M_PI * x[1]);
109+
Real sin_pi_z_2 = std::sin(M_PI * x[2]);
110+
sin_pi_x_2 *= sin_pi_x_2;
111+
sin_pi_y_2 *= sin_pi_y_2;
112+
sin_pi_z_2 *= sin_pi_z_2;
113+
114+
Real sin_2_pi_x = std::sin(2.0 * M_PI * x[0]);
115+
Real sin_2_pi_y = std::sin(2.0 * M_PI * x[1]);
116+
Real sin_2_pi_z = std::sin(2.0 * M_PI * x[2]);
117+
return {{
118+
2.0 * sin_pi_x_2 * sin_2_pi_y * sin_2_pi_z,
119+
-sin_2_pi_x * sin_pi_y_2 * sin_2_pi_z,
120+
-sin_2_pi_x * sin_2_pi_y * sin_pi_z_2}};
121+
122+
}
123+
124+
void
125+
advect()
126+
{
127+
#ifdef _MSC_VER
128+
concurrency::parallel_for_each
129+
#else
130+
__gnu_parallel::for_each
131+
#endif
132+
(positions.begin(), positions.end(), [&](std::array<Real, 3>& x)
133+
{
134+
std::array<Real, 3> v = enright_velocity_field(x);
135+
x[0] += 0.005 * v[0];
136+
x[1] += 0.005 * v[1];
137+
x[2] += 0.005 * v[1];
138+
}
139+
);
140+
}
141+
142+
int main(int argc, char* argv[])
143+
{
144+
Real min_x = std::numeric_limits<Real>::max();
145+
Real max_x = std::numeric_limits<Real>::min();
146+
positions.reserve(N * N * N);
147+
for (unsigned int i = 0; i < N; ++i)
148+
{
149+
for (unsigned int j = 0; j < N; ++j)
150+
{
151+
for (unsigned int k = 0; k < N; ++k)
152+
{
153+
std::array<Real, 3> x = {{
154+
r_omega * (2.0 * static_cast<Real>(i) / static_cast<Real>(N-1)-1.0),
155+
r_omega * (2.0 * static_cast<Real>(j) / static_cast<Real>(N-1)-1.0),
156+
r_omega * (2.0 * static_cast<Real>(k) / static_cast<Real>(N-1)-1.0)}};
157+
158+
Real l2 = x[0] * x[0] + x[1] * x[1] + x[2] * x[2];
159+
if (l2 < r_omega2)
160+
{
161+
x[0] += 0.35;
162+
x[1] += 0.35;
163+
x[2] += 0.35;
164+
positions.push_back(x);
165+
if (min_x > x[0])
166+
{
167+
min_x = x[0];
168+
}
169+
if (max_x < x[0])
170+
{
171+
max_x = x[0];
172+
}
173+
}
174+
}
175+
}
176+
}
177+
std::random_shuffle(positions.begin(), positions.end());
178+
179+
NeighborhoodSearch nsearch(radius, true);
180+
nsearch.add_point_set(positions.front().data(), positions.size(), true, true);
181+
nsearch.find_neighbors();
182+
183+
std::cout << "#Points = " << positions.size() << std::endl;
184+
std::cout << "Search radius = " << radius << std::endl;
185+
std::cout << "Min x = " << min_x << std::endl;
186+
std::cout << "Max x = " << max_x << std::endl;
187+
std::cout << "Average number of neighbors = " << compute_average_number_of_neighbors(nsearch) << std::endl;
188+
std::cout << "Average index distance prior to z-sort = " << compute_average_distance(nsearch) << std::endl;
189+
190+
nsearch.z_sort();
191+
for (auto const& d : nsearch.point_sets())
192+
{
193+
d.sort_field(positions.data());
194+
}
195+
nsearch.find_neighbors();
196+
//compare_with_bruteforce_search(nsearch);
197+
198+
std::cout << "Average index distance after z-sort = " << compute_average_distance(nsearch) << std::endl;
199+
200+
std::cout << "Moving points:" << std::endl;
201+
for (int i = 0; i < N_enright_steps; ++i)
202+
{
203+
std::cout << "Enright step " << i << ". ";
204+
advect();
205+
auto t0 = std::chrono::high_resolution_clock::now();
206+
nsearch.find_neighbors();
207+
std::cout << "Neighborhood search took " << std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::high_resolution_clock::now() - t0).count() << "ms" << std::endl;
208+
//compare_with_bruteforce_search(nsearch);
209+
}
210+
211+
return 0;
212+
}

‎extern/libmorton/LICENSE

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
MIT License
2+
3+
Copyright (c) 2016 Jeroen Baert
4+
5+
Permission is hereby granted, free of charge, to any person obtaining a copy
6+
of this software and associated documentation files (the "Software"), to deal
7+
in the Software without restriction, including without limitation the rights
8+
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
9+
copies of the Software, and to permit persons to whom the Software is
10+
furnished to do so, subject to the following conditions:
11+
12+
The above copyright notice and this permission notice shall be included in all
13+
copies or substantial portions of the Software.
14+
15+
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16+
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17+
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18+
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19+
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20+
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
21+
SOFTWARE.

0 commit comments

Comments
 (0)