Fix cmake gmsh (#110)

* Update FindCmake

.circleci/config.yml

@@ -42,6 +42,27 @@ jobs:
             - store_artifacts:
                 path: final
                 destination: packages
+    build_pcl:
+        docker:
+          - image: diehlpk/nonlocalmodels:baseimage
+        steps:
+            - checkout
+            - run:
+                name: Configure 
+                command: | 
+                  mkdir -p build 
+                  cd build 
+                  cmake -DCMAKE_BUILD_TYPE=Release -DEnable_Documentation=ON -DEnable_Tools=ON -DEnable_PCL=ON ..
+            - run:
+                name: Build
+                command: |  
+                  cd build 
+                  make -j 2
+            - run:
+                name: Test
+                command: |
+                  cd build 
+                  make test
     codecoverage:
       docker: 
         - image: diehlpk/nonlocalmodels:baseimage
@@ -155,4 +176,9 @@ workflows:
                 ignore: gh-pages
             requires:
               - build
-
+      - build_pcl:
+            filters:
+              branches:
+                ignore: gh-pages
+            requires:
+              - build

CMakeLists.txt

@@ -95,6 +95,12 @@ find_package(BlazeIterative REQUIRED)
 find_package(LAPACK REQUIRED)
 find_package(BLAS REQUIRED)
 
+
+set(Enable_PCL FALSE CACHE BOOL "Enables PCL for faster neighborsearch")
+if(${Enable_PCL})
+    find_package(PCL REQUIRED)
+endif()
+
 #####################################################################
 # Set the compiler flags
 #####################################################################

README.md

@@ -23,9 +23,13 @@ as the build system and has following dependencies:
   * [Blaze_Iterative](https://github.com/STEllAR-GROUP/BlazeIterative) master, 
   * [Boost](https://www.boost.org/) 1.73,
   * [gmsh](https://gmsh.info/) 4.6,
-  * [VTK](https://www.vtk.org) 8.2, and 
+  * [VTK](https://www.vtk.org) 9.0, and 
   * [YAML-CPP](https://github.com/jbeder/yaml-cpp) 0.6.3.
 
+Following dependencies are optional, but are recommend for large simualtions:
+
+  * [PCL](https://github.com/PointCloudLibrary/pcl) 1.11. 
+
 Note that the above provided versions are the minimal rquired versions and we recommend to use these to build the code. 
 
 We provide following support to build the code

cmake/FindGmsh.cmake

@@ -6,17 +6,17 @@
 find_package(PkgConfig)
 
 find_library(GMSH_LIB
-        NAMES gmsh
-        PATH /usr/lib64 /usr/local/lib64 /usr/lib/ /usr/local/lib "${GMSH_DIR}/lib/"
-	PATH_SUFFIXES lib lib64)
+        NAMES libgmsh.a libgmsh.so
+        HINTS /usr/lib64 /usr/local/lib64 /usr/lib/ /usr/local/lib "${GMSH_DIR}"
+	    PATH_SUFFIXES lib lib64)
 find_path(GMSH_INCLUDE gmsh.h HINTS /usr/include /usr/local/include "${GMSH_DIR}/include/")
 
 mark_as_advanced(GMSH_DIR)
 mark_as_advanced(GMSH_LIB)
 mark_as_advanced(GMSH_INCLUDE)
 
 if (NOT GMSH_LIB)
-    message(FATAL_ERROR "Gmsh Library not found: Specify the GMSH_DIR where yaml cpp is located")
+    message(FATAL_ERROR "Gmsh Library not found: Specify the GMSH_DIR where Gmsh is located")
 else ()
     include_directories(${GMSH_INCLUDE})
 endif ()

docs/content/cmake-options.md

@@ -6,6 +6,7 @@
 * Enable_Tools : Enables the tools to the build target (Default = False)
 * Enable_Expensive_Tests : Enables the computation intense tests (Default = False)
 * Enable_RPM: Enables to generate RPM packages (Default = False)
+* Enable_PCL: Add [PCL](https://github.com/PointCloudLibrary/pcl) for faster neighborseach. Recommend for large amount of nodes (Default=False)
 
 ## General options
 

src/geometry/CMakeLists.txt

@@ -6,5 +6,15 @@
 
 AUX_SOURCE_DIRECTORY(./ SOURCES)
 
+if(${Enable_PCL})
+    add_definitions(-DENABLE_PCL)
+    include_directories(${PCL_INCLUDE_DIRS})
+    link_directories(${PCL_LIBRARY_DIRS})
+endif()
+
 add_hpx_library(Geometry
     SOURCES ${SOURCES})
+
+if(${Enable_PCL})
+    target_link_libraries(Geometry PUBLIC ${PCL_LIBRARIES})
+endif()

src/geometry/neighbor.cpp

@@ -14,11 +14,58 @@
 #include "util/compare.h"
 #include "util/utilIO.h"
 
+#ifdef ENABLE_PCL
+
+#include <pcl/point_cloud.h>
+#include <pcl/kdtree/kdtree_flann.h>
+
+#endif
+
 geometry::Neighbor::Neighbor(const double &horizon, inp::NeighborDeck *deck,
                              const std::vector<util::Point3> *nodes)
     : d_neighborDeck_p(deck) {
   d_neighbors.resize(nodes->size());
 
+#ifdef ENABLE_PCL
+
+  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>);
+
+  cloud->points.resize(nodes->size());
+
+  hpx::for_loop(hpx::parallel::execution::par, 0, nodes->size(),
+                [this, cloud, nodes](boost::uint64_t i) {
+                  (*cloud)[i].x = (*nodes)[i].d_x;
+                  (*cloud)[i].y = (*nodes)[i].d_y;
+                  (*cloud)[i].z = (*nodes)[i].d_z;
+                });
+
+  pcl::KdTreeFLANN<pcl::PointXYZ> kdtree;
+
+  kdtree.setInputCloud(cloud);
+
+  hpx::for_loop(hpx::parallel::execution::par, 0, nodes->size(),
+                [this, kdtree, nodes, horizon](boost::uint64_t i) {
+                  std::vector<int> neighs;
+                  std::vector<float> pointRadiusSquaredDistance;
+
+                  pcl::PointXYZ searchPoint;
+                  searchPoint.x = (*nodes)[i].d_x;
+                  searchPoint.y = (*nodes)[i].d_y;
+                  searchPoint.z = (*nodes)[i].d_z;
+
+                  this->d_neighbors[i] = std::vector<size_t>();
+
+                  if (kdtree.radiusSearch(searchPoint, horizon, neighs,
+                                          pointRadiusSquaredDistance) > 0) {
+                    for (std::size_t j = 0; j < neighs.size(); ++j)
+                      if (neighs[j] != i) {
+                        this->d_neighbors[i].push_back(size_t(neighs[j]));
+                      }
+                  }
+                });
+
+#else
+
   auto f = hpx::for_loop(
       hpx::parallel::execution::par(hpx::parallel::execution::task), 0,
       nodes->size(), [this, horizon, nodes](boost::uint64_t i) {
@@ -40,6 +87,8 @@ geometry::Neighbor::Neighbor(const double &horizon, inp::NeighborDeck *deck,
       });  // end of parallel for loop
 
   f.get();
+
+#endif
 }
 
 const std::vector<size_t> &geometry::Neighbor::getNeighbors(const size_t &i) {

src/rw/vtkWriter.cpp

@@ -70,7 +70,7 @@ void rw::writer::VtkWriter::appendMesh(const std::vector<util::Point3> *nodes,
   // element type
   int *cell_types = new int[num_elems];
 
-  vtkIdType ids[num_vertex];
+  vtkIdType *ids = new vtkIdType[num_vertex];
 
   for (size_t i = 0; i < num_elems; i++) {
     // get ids of vertex of this element