Backup before updating libtimestep

egor-demidov · Mar 4, 2024 · f850dc9 · f850dc9
1 parent 2e2508a
commit f850dc9
Show file tree

Hide file tree

Showing 7 changed files with 218 additions and 65 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -59,6 +59,8 @@ add_executable(sintered_test test/sintered.cpp test/mass_distribution.cpp)
 add_executable(alt_sintered_test test/alt_sintered.cpp writer.cpp test/mass_distribution.cpp)
 # Particle is collided with a plane
 add_executable(surface_test test/surface.cpp writer.cpp test/compute_energy.cpp)
+# Same as hamaker_test, but with truncation and neighbor tracking
+add_executable(hamaker_neighbor_test test/hamaker_neighbor.cpp)
 
 add_executable(sintered_2_test test/sintered_2.cpp writer.cpp)
 

diff --git a/include/libgran/granular_system/binary_force_container.h b/include/libgran/granular_system/binary_force_container.h
@@ -0,0 +1,38 @@
+//
+// Created by egor on 3/4/24.
+//
+
+#ifndef LIBGRAN_BINARY_FORCE_CONTAINER_H
+#define LIBGRAN_BINARY_FORCE_CONTAINER_H
+
+#include <vector>
+
+template <
+        typename field_value_t,
+        typename real_t,
+        typename... binary_force_functors_t>
+struct binary_force_functor_container {
+    typedef std::vector<field_value_t> field_container_t;
+
+    binary_force_functor_container(binary_force_functors_t & ... force_functors)
+            : binary_force_functors{force_functors...} {}
+
+    std::pair<field_value_t, field_value_t> operator () (size_t i, size_t j, field_container_t const & x,
+                                                         field_container_t const & v,
+                                                         field_container_t const & theta,
+                                                         field_container_t const & omega,
+                                                         real_t t) {
+        static_assert(std::tuple_size<decltype(binary_force_functors)>::value > 0, "at least one binary force functor must be provided "
+                                                                                   "as a template parameter");
+
+        std::pair<field_value_t, field_value_t> accelerations = std::apply([i, j, &x, &v, &theta, &omega, t] (auto & ... e) -> std::pair<field_value_t, field_value_t> {
+            return (e(i, j, x, v, theta, omega, t) + ...);
+        }, binary_force_functors);
+        return accelerations;
+    }
+
+    std::tuple<binary_force_functors_t & ...> binary_force_functors;
+};
+
+
+#endif //LIBGRAN_BINARY_FORCE_CONTAINER_H
diff --git a/include/libgran/granular_system/granular_system.h b/include/libgran/granular_system/granular_system.h
@@ -18,71 +18,9 @@
 #pragma message("libgran: using OpenMP parallelization")
 #endif //LIBGRAN_USE_OMP
 
-template<typename T1, typename T2>
-std::pair<T1, T2> operator + (std::pair<T1, T2> const & lhs, std::pair<T1, T2> const & rhs) {
-    return std::make_pair(lhs.first + rhs.first, lhs.second + rhs.second);
-}
-
-template<typename T1, typename T2>
-std::pair<T1, T2> & operator += (std::pair<T1, T2> & lhs, std::pair<T1, T2> const & rhs) {
-    lhs.first += rhs.first;
-    lhs.second += rhs.second;
-    return lhs;
-}
-
-template <
-    typename field_value_t,
-    typename real_t,
-    typename... binary_force_functors_t>
-struct binary_force_functor_container {
-    typedef std::vector<field_value_t> field_container_t;
-
-    binary_force_functor_container(binary_force_functors_t & ... force_functors)
-            : binary_force_functors{force_functors...} {}
-
-    std::pair<field_value_t, field_value_t> operator () (size_t i, size_t j, field_container_t const & x,
-                                                                  field_container_t const & v,
-                                                                  field_container_t const & theta,
-                                                                  field_container_t const & omega,
-                                                                  real_t t) {
-        static_assert(std::tuple_size<decltype(binary_force_functors)>::value > 0, "at least one binary force functor must be provided "
-                                                                            "as a template parameter");
-
-        std::pair<field_value_t, field_value_t> accelerations = std::apply([i, j, &x, &v, &theta, &omega, t] (auto & ... e) -> std::pair<field_value_t, field_value_t> {
-            return (e(i, j, x, v, theta, omega, t) + ...);
-        }, binary_force_functors);
-        return accelerations;
-    }
-
-    std::tuple<binary_force_functors_t & ...> binary_force_functors;
-};
-
-template <
-    typename field_value_t,
-    typename real_t,
-    typename... unary_force_functors_t>
-struct unary_force_functor_container {
-    typedef std::vector<field_value_t> field_container_t;
-
-    unary_force_functor_container(unary_force_functors_t & ... force_functors)
-            : unary_force_functors{force_functors...} {}
-
-    std::pair<field_value_t, field_value_t> operator () (size_t i, field_container_t const & x,
-                                                                  field_container_t const & v,
-                                                                  field_container_t const & theta,
-                                                                  field_container_t const & omega,
-                                                                  real_t t) {
-        static_assert(std::tuple_size<decltype(unary_force_functors)>::value > 0, "at least one unary force functor must be provided "
-                                                                            "as a template parameter");
-
-        std::pair<field_value_t, field_value_t> accelerations = std::apply([i, &x, &v, &theta, &omega, t] (auto & ... e) -> std::pair<field_value_t, field_value_t> {
-            return (e(i, x, v, theta, omega, t) + ...);
-        }, unary_force_functors);
-        return accelerations;
-    }
-
-    std::tuple<unary_force_functors_t & ...> unary_force_functors;
-};
+#include "../utility/operator_overloads.h"
+#include "binary_force_container.h"
+#include "unary_force_container.h"
 
 template <
     typename field_value_t,

diff --git a/include/libgran/granular_system/granular_system_neighbor_list.h b/include/libgran/granular_system/granular_system_neighbor_list.h
@@ -0,0 +1,92 @@
+//
+// Created by egor on 3/4/24.
+//
+
+#ifndef LIBGRAN_GRANULAR_SYSTEM_NEIGHBOR_LIST_H
+#define LIBGRAN_GRANULAR_SYSTEM_NEIGHBOR_LIST_H
+
+
+#include <libtimestep/rotational_integrator/rotational_integrator.h>
+#include <libtimestep/rotational_step_handler/rotational_step_handler.h>
+
+#ifndef LIBGRAN_USE_OMP
+#define binary_system_implementation rotational_binary_system
+#include <libtimestep/rotational_system/rotational_binary_system.h>
+#pragma message("libgran: using C++ 17 parallel algorithms")
+#else
+#define binary_system_implementation rotational_binary_system_omp
+#include <libtimestep/rotational_system/rotational_binary_system_omp.h>
+#pragma message("libgran: using OpenMP parallelization")
+#endif //LIBGRAN_USE_OMP
+
+#include "../utility/operator_overloads.h"
+#include "binary_force_container.h"
+#include "unary_force_container.h"
+
+template <
+        typename field_value_t,
+        typename real_t,
+        template <
+        typename _field_container_t,
+        typename _field_value_t,
+        typename _real_t,
+        typename _functor_t,
+        template <
+        typename __field_container_t,
+        typename __field_value_t>
+        typename _step_handler_t>
+        typename integrator_t,
+        template <
+        typename _field_container_t,
+        typename _field_value_t>
+        typename step_handler_t,
+        typename binary_force_functor_container_t,
+        typename unary_force_functor_container_t>
+class granular_system_neignbor_list : public binary_system_implementation<field_value_t, real_t, integrator_t, step_handler_t,
+        granular_system_neignbor_list<field_value_t, real_t, integrator_t, step_handler_t,
+                binary_force_functor_container_t, unary_force_functor_container_t>,
+        (std::tuple_size<decltype(unary_force_functor_container_t::unary_force_functors)>::value > 0)> {
+public:
+    typedef std::vector<field_value_t> field_container_t;
+
+    granular_system_neignbor_list(granular_system_neignbor_list const &) = delete;
+
+    granular_system_neignbor_list(field_container_t x0, field_container_t v0,
+                    field_container_t theta0, field_container_t omega0,
+                    real_t t0, field_value_t field_zero, real_t real_zero, step_handler_t<field_container_t, field_value_t> & step_handler,
+                    binary_force_functor_container_t binary_force_functors, unary_force_functor_container_t unary_force_functors) :
+            binary_system_implementation<field_value_t, real_t, integrator_t, step_handler_t,
+                    granular_system_neignbor_list<field_value_t, real_t, integrator_t, step_handler_t, binary_force_functor_container_t, unary_force_functor_container_t>,
+                    (std::tuple_size<decltype(unary_force_functor_container_t::unary_force_functors)>::value > 0)>
+                    (std::move(x0), std::move(v0), std::move(theta0),
+                     std::move(omega0), t0, field_zero, real_zero, *this, step_handler),
+            binary_force_functors(std::move(binary_force_functors)),
+            unary_force_functors(std::move(unary_force_functors)) {}
+
+    std::pair<field_value_t, field_value_t> compute_accelerations(size_t i, size_t j, field_container_t const & x,
+                                                                  field_container_t const & v,
+                                                                  field_container_t const & theta,
+                                                                  field_container_t const & omega,
+                                                                  real_t t) {
+
+        return binary_force_functors(i, j, x, v, theta, omega, t);
+    }
+
+    std::pair<field_value_t, field_value_t> compute_accelerations(size_t i, field_container_t const & x,
+                                                                  field_container_t const & v,
+                                                                  field_container_t const & theta,
+                                                                  field_container_t const & omega,
+                                                                  real_t t) {
+
+        return unary_force_functors(i, x, v, theta, omega, t);
+    }
+
+private:
+    const real_t r_verlet;  // Verlet radius for neighbor list computation
+
+    binary_force_functor_container_t binary_force_functors;
+    unary_force_functor_container_t unary_force_functors;
+};
+
+
+#endif //LIBGRAN_GRANULAR_SYSTEM_NEIGHBOR_LIST_H
diff --git a/include/libgran/granular_system/unary_force_container.h b/include/libgran/granular_system/unary_force_container.h
@@ -0,0 +1,37 @@
+//
+// Created by egor on 3/4/24.
+//
+
+#ifndef LIBGRAN_UNARY_FORCE_CONTAINER_H
+#define LIBGRAN_UNARY_FORCE_CONTAINER_H
+
+#include <vector>
+
+template <
+        typename field_value_t,
+        typename real_t,
+        typename... unary_force_functors_t>
+struct unary_force_functor_container {
+    typedef std::vector<field_value_t> field_container_t;
+
+    unary_force_functor_container(unary_force_functors_t & ... force_functors)
+            : unary_force_functors{force_functors...} {}
+
+    std::pair<field_value_t, field_value_t> operator () (size_t i, field_container_t const & x,
+                                                         field_container_t const & v,
+                                                         field_container_t const & theta,
+                                                         field_container_t const & omega,
+                                                         real_t t) {
+        static_assert(std::tuple_size<decltype(unary_force_functors)>::value > 0, "at least one unary force functor must be provided "
+                                                                                  "as a template parameter");
+
+        std::pair<field_value_t, field_value_t> accelerations = std::apply([i, &x, &v, &theta, &omega, t] (auto & ... e) -> std::pair<field_value_t, field_value_t> {
+            return (e(i, x, v, theta, omega, t) + ...);
+        }, unary_force_functors);
+        return accelerations;
+    }
+
+    std::tuple<unary_force_functors_t & ...> unary_force_functors;
+};
+
+#endif //LIBGRAN_UNARY_FORCE_CONTAINER_H
diff --git a/include/libgran/utility/operator_overloads.h b/include/libgran/utility/operator_overloads.h
@@ -0,0 +1,22 @@
+//
+// Created by egor on 3/4/24.
+//
+
+#ifndef LIBGRAN_OPERATOR_OVERLOADS_H
+#define LIBGRAN_OPERATOR_OVERLOADS_H
+
+#include <utility>
+
+template<typename T1, typename T2>
+std::pair<T1, T2> operator + (std::pair<T1, T2> const & lhs, std::pair<T1, T2> const & rhs) {
+return std::make_pair(lhs.first + rhs.first, lhs.second + rhs.second);
+}
+
+template<typename T1, typename T2>
+std::pair<T1, T2> & operator += (std::pair<T1, T2> & lhs, std::pair<T1, T2> const & rhs) {
+    lhs.first += rhs.first;
+    lhs.second += rhs.second;
+    return lhs;
+}
+
+#endif //LIBGRAN_OPERATOR_OVERLOADS_H
diff --git a/test/hamaker_neighbor.cpp b/test/hamaker_neighbor.cpp
@@ -0,0 +1,24 @@
+//
+// Created by egor on 3/4/24.
+//
+
+#include <iostream>
+#include <vector>
+#include <chrono>
+
+#ifdef _GNU_SOURCE
+#include <cfenv>
+#define enable_fp_exceptions() feenableexcept(FE_INVALID | FE_OVERFLOW | FE_DIVBYZERO)
+#elif defined(_MSC_VER)
+#define enable_fp_exceptions()
+#elif defined(__APPLE__)
+#define enable_fp_exceptions()
+#else
+#error "Unsupported system"
+#endif
+
+#include <Eigen/Eigen>
+
+#include <libgran/contact_force/contact_force.h>
+#include <libgran/hamaker_force/hamaker_force.h>
+#include <libgran/granular_system/granular_system_neighbor_list.h>