Adapted to the changes in libtimestep

include/libgran/granular_system/granular_system.h

@@ -14,6 +14,60 @@ std::pair<T1, T2> operator + (std::pair<T1, T2> const & lhs, std::pair<T1, T2> c
     return std::make_pair(lhs.first + rhs.first, lhs.second + rhs.second);
 }
 
+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 binary 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;
+};
+
 template <
     typename field_value_t,
     typename real_t,
@@ -31,9 +85,11 @@ template <
         typename _field_container_t,
         typename _field_value_t>
     typename step_handler_t,
-    typename... force_functors_t>
+    typename binary_force_functor_container_t,
+    typename unary_force_functor_container_t>
 class granular_system : public rotational_binary_system<field_value_t, real_t, integrator_t, step_handler_t,
-        granular_system<field_value_t, real_t, integrator_t, step_handler_t, force_functors_t...>> {
+        granular_system<field_value_t, real_t, integrator_t, step_handler_t,
+        binary_force_functor_container_t, unary_force_functor_container_t>> {
 public:
     typedef std::vector<field_value_t> field_container_t;
 
@@ -42,28 +98,35 @@ class granular_system : public rotational_binary_system<field_value_t, real_t, i
     granular_system(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,
-                    force_functors_t & ... force_functors) :
+                    binary_force_functor_container_t binary_force_functors, unary_force_functor_container_t unary_force_functors) :
             rotational_binary_system<field_value_t, real_t, integrator_t, step_handler_t,
-                    granular_system<field_value_t, real_t, integrator_t, step_handler_t, force_functors_t...>>
+                    granular_system<field_value_t, real_t, integrator_t, step_handler_t, binary_force_functor_container_t, unary_force_functor_container_t>>
                     (std::move(x0), std::move(v0), std::move(theta0),
-                     std::move(omega0), t0, field_zero, real_zero, *this, step_handler), force_functors{force_functors...} {}
+                     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) {
-        static_assert(std::tuple_size<decltype(force_functors)>::value > 0, "at least one 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) + ...);
-        }, force_functors);
-        return accelerations;
+        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:
-    std::tuple<force_functors_t & ...> force_functors;
+    binary_force_functor_container_t binary_force_functors;
+    unary_force_functor_container_t unary_force_functors;
 };
 
 #endif //LIBGRAN_GRANULAR_SYSTEM_H

include/libgran/no_unary_force/no_unary_force.h

@@ -0,0 +1,26 @@
+//
+// Created by egor on 1/28/24.
+//
+
+#ifndef LIBGRAN_NO_UNARY_FORCE_H
+#define LIBGRAN_NO_UNARY_FORCE_H
+
+template <typename field_value_t, typename real_t>
+struct no_unary_force_functor {
+    explicit no_unary_force_functor(field_value_t zero_field) : zero_field(std::move(zero_field)) {}
+
+    std::pair<field_value_t, field_value_t> operator () (size_t i [[maybe_unused]],
+                                                         std::vector<field_value_t> const & x [[maybe_unused]],
+                                                         std::vector<field_value_t> const & v [[maybe_unused]],
+                                                         std::vector<field_value_t> const & theta [[maybe_unused]],
+                                                         std::vector<field_value_t> const & omega [[maybe_unused]],
+                                                         real_t t [[maybe_unused]]) const {
+
+        return std::make_pair(zero_field, zero_field);
+    }
+
+private:
+    const field_value_t zero_field;
+};
+
+#endif //LIBGRAN_NO_UNARY_FORCE_H

main.cpp

@@ -10,6 +10,7 @@
 
 #include <libgran/contact_force/contact_force.h>
 #include <libgran/hamaker_force/hamaker_force.h>
+#include <libgran/no_unary_force/no_unary_force.h>
 #include <libgran/granular_system/granular_system.h>
 
 #include "writer.h"
@@ -100,14 +101,29 @@ int main() {
     hamaker_functor<Eigen::Vector3d, double> hamaker_model(A, h0,
                                                            r_part, mass, Eigen::Vector3d::Zero(), 0.0);
 
+    binary_force_functor_container<Eigen::Vector3d, double,
+        contact_force_functor<Eigen::Vector3d, double>,
+        hamaker_functor<Eigen::Vector3d, double>>
+            binary_force_functors{contact_force_model, hamaker_model};
+
+    no_unary_force_functor<Eigen::Vector3d, double> no_unary_force(Eigen::Vector3d::Zero());
+
+    unary_force_functor_container<Eigen::Vector3d, double,
+        no_unary_force_functor<Eigen::Vector3d, double>>
+            unary_force_functors(no_unary_force);
+
     // Create an instance of granular_system using the contact force model
     // Using velocity Verlet integrator for rotational systems and a default
     // step handler for rotational systems
     granular_system<Eigen::Vector3d, double, rotational_velocity_verlet_half,
-            rotational_step_handler, contact_force_functor<Eigen::Vector3d, double>,
-            hamaker_functor<Eigen::Vector3d, double>> system(x0,
-                                                             v0, theta0, omega0, 0.0, Eigen::Vector3d::Zero(), 0.0,
-                                                             step_handler_instance, contact_force_model, hamaker_model);
+        rotational_step_handler,
+        binary_force_functor_container<Eigen::Vector3d, double,
+                contact_force_functor<Eigen::Vector3d, double>,
+                hamaker_functor<Eigen::Vector3d, double>>,
+        unary_force_functor_container<Eigen::Vector3d, double,
+                no_unary_force_functor<Eigen::Vector3d, double>>> system(x0,
+            v0, theta0, omega0, 0.0, Eigen::Vector3d::Zero(), 0.0,
+            step_handler_instance, binary_force_functors, unary_force_functors);
 
 
     // Output stream for data

test/contact.cpp

@@ -9,6 +9,7 @@
 #include <Eigen/Eigen>
 
 #include <libgran/contact_force/contact_force.h>
+#include <libgran/no_unary_force/no_unary_force.h>
 #include <libgran/granular_system/granular_system.h>
 
 #include "compute_energy.h"
@@ -73,12 +74,21 @@ int main() {
        k, gamma_n, k, gamma_t, mu, phi, k, gamma_r, mu_o, phi, k, gamma_o, mu_o, phi,
        r_part, mass, inertia, dt, Eigen::Vector3d::Zero(), 0.0);
 
+    binary_force_functor_container<Eigen::Vector3d, double, contact_force_functor<Eigen::Vector3d, double>>
+        binary_force_functors{contact_force_model};
+
+    no_unary_force_functor<Eigen::Vector3d, double> no_unary_force(Eigen::Vector3d::Zero());
+
+    unary_force_functor_container<Eigen::Vector3d, double, no_unary_force_functor<Eigen::Vector3d, double>>
+        unary_force_functors(no_unary_force);
+
     // Create an instance of granular_system using the contact force model
     // Using velocity Verlet integrator for rotational systems and a default
     // step handler for rotational systems
     granular_system<Eigen::Vector3d, double, rotational_velocity_verlet_half,
-        rotational_step_handler, contact_force_functor<Eigen::Vector3d, double>> system(x0,
-            v0, theta0, omega0, 0.0, Eigen::Vector3d::Zero(), 0.0, step_handler_instance, contact_force_model);
+        rotational_step_handler, binary_force_functor_container<Eigen::Vector3d, double, contact_force_functor<Eigen::Vector3d, double>>,
+            unary_force_functor_container<Eigen::Vector3d, double, no_unary_force_functor<Eigen::Vector3d, double>>> system(x0,
+            v0, theta0, omega0, 0.0, Eigen::Vector3d::Zero(), 0.0, step_handler_instance, binary_force_functors, unary_force_functors);
 
     auto start_time = std::chrono::high_resolution_clock::now();
     for (size_t n = 0; n < n_steps; n ++) {

test/hamaker.cpp

@@ -10,6 +10,7 @@
 
 #include <libgran/contact_force/contact_force.h>
 #include <libgran/hamaker_force/hamaker_force.h>
+#include <libgran/no_unary_force/no_unary_force.h>
 #include <libgran/granular_system/granular_system.h>
 
 #include "compute_energy.h"
@@ -83,14 +84,29 @@ int main() {
     hamaker_functor<Eigen::Vector3d, double> hamaker_model(A, h0,
        r_part, mass, Eigen::Vector3d::Zero(), 0.0);
 
+    binary_force_functor_container<Eigen::Vector3d, double,
+        contact_force_functor<Eigen::Vector3d, double>,
+        hamaker_functor<Eigen::Vector3d, double>>
+            binary_force_functors{contact_force_model, hamaker_model};
+
+    no_unary_force_functor<Eigen::Vector3d, double> no_unary_force(Eigen::Vector3d::Zero());
+
+    unary_force_functor_container<Eigen::Vector3d, double,
+        no_unary_force_functor<Eigen::Vector3d, double>>
+            unary_force_functors(no_unary_force);
+
     // Create an instance of granular_system using the contact force model
     // Using velocity Verlet integrator for rotational systems and a default
     // step handler for rotational systems
     granular_system<Eigen::Vector3d, double, rotational_velocity_verlet_half,
-        rotational_step_handler, contact_force_functor<Eigen::Vector3d, double>,
-        hamaker_functor<Eigen::Vector3d, double>> system(x0,
+        rotational_step_handler,
+        binary_force_functor_container<Eigen::Vector3d, double,
+                contact_force_functor<Eigen::Vector3d, double>,
+                hamaker_functor<Eigen::Vector3d, double>>,
+        unary_force_functor_container<Eigen::Vector3d, double,
+                no_unary_force_functor<Eigen::Vector3d, double>>> system(x0,
             v0, theta0, omega0, 0.0, Eigen::Vector3d::Zero(), 0.0,
-            step_handler_instance, contact_force_model, hamaker_model);
+            step_handler_instance, binary_force_functors, unary_force_functors);
 
     auto start_time = std::chrono::high_resolution_clock::now();
     for (size_t n = 0; n < n_steps; n ++) {

test/hamaker_2.cpp

@@ -10,6 +10,7 @@
 
 #include <libgran/contact_force/contact_force.h>
 #include <libgran/hamaker_force/hamaker_force.h>
+#include <libgran/no_unary_force/no_unary_force.h>
 #include <libgran/granular_system/granular_system.h>
 
 #include "../writer.h"
@@ -91,14 +92,29 @@ int main() {
     hamaker_functor<Eigen::Vector3d, double> hamaker_model(A, h0,
        r_part, mass, Eigen::Vector3d::Zero(), 0.0);
 
+    binary_force_functor_container<Eigen::Vector3d, double,
+        contact_force_functor<Eigen::Vector3d, double>,
+        hamaker_functor<Eigen::Vector3d, double>>
+            binary_force_functors{contact_force_model, hamaker_model};
+
+    no_unary_force_functor<Eigen::Vector3d, double> no_unary_force(Eigen::Vector3d::Zero());
+
+    unary_force_functor_container<Eigen::Vector3d, double,
+        no_unary_force_functor<Eigen::Vector3d, double>>
+            unary_force_functors(no_unary_force);
+
     // Create an instance of granular_system using the contact force model
     // Using velocity Verlet integrator for rotational systems and a default
     // step handler for rotational systems
     granular_system<Eigen::Vector3d, double, rotational_velocity_verlet_half,
-        rotational_step_handler, contact_force_functor<Eigen::Vector3d, double>,
-        hamaker_functor<Eigen::Vector3d, double>> system(x0,
+        rotational_step_handler,
+        binary_force_functor_container<Eigen::Vector3d, double,
+                contact_force_functor<Eigen::Vector3d, double>,
+                hamaker_functor<Eigen::Vector3d, double>>,
+        unary_force_functor_container<Eigen::Vector3d, double,
+                no_unary_force_functor<Eigen::Vector3d, double>>> system(x0,
             v0, theta0, omega0, 0.0, Eigen::Vector3d::Zero(), 0.0,
-            step_handler_instance, contact_force_model, hamaker_model);
+            step_handler_instance, binary_force_functors, unary_force_functors);
 
     auto start_time = std::chrono::high_resolution_clock::now();
     for (size_t n = 0; n < n_steps; n ++) {

test/hamaker_3.cpp

@@ -10,6 +10,7 @@
 
 #include <libgran/contact_force/contact_force.h>
 #include <libgran/hamaker_force/hamaker_force.h>
+#include <libgran/no_unary_force/no_unary_force.h>
 #include <libgran/granular_system/granular_system.h>
 
 #include "../writer.h"
@@ -97,14 +98,29 @@ int main() {
     hamaker_functor<Eigen::Vector3d, double> hamaker_model(A, h0,
        r_part, mass, Eigen::Vector3d::Zero(), 0.0);
 
+    binary_force_functor_container<Eigen::Vector3d, double,
+        contact_force_functor<Eigen::Vector3d, double>,
+        hamaker_functor<Eigen::Vector3d, double>>
+            binary_force_functors{contact_force_model, hamaker_model};
+
+    no_unary_force_functor<Eigen::Vector3d, double> no_unary_force(Eigen::Vector3d::Zero());
+
+    unary_force_functor_container<Eigen::Vector3d, double,
+        no_unary_force_functor<Eigen::Vector3d, double>>
+            unary_force_functors(no_unary_force);
+
     // Create an instance of granular_system using the contact force model
     // Using velocity Verlet integrator for rotational systems and a default
     // step handler for rotational systems
     granular_system<Eigen::Vector3d, double, rotational_velocity_verlet_half,
-        rotational_step_handler, contact_force_functor<Eigen::Vector3d, double>,
-        hamaker_functor<Eigen::Vector3d, double>> system(x0,
+        rotational_step_handler,
+        binary_force_functor_container<Eigen::Vector3d, double,
+                contact_force_functor<Eigen::Vector3d, double>,
+                hamaker_functor<Eigen::Vector3d, double>>,
+        unary_force_functor_container<Eigen::Vector3d, double,
+                no_unary_force_functor<Eigen::Vector3d, double>>> system(x0,
             v0, theta0, omega0, 0.0, Eigen::Vector3d::Zero(), 0.0,
-            step_handler_instance, contact_force_model, hamaker_model);
+            step_handler_instance, binary_force_functors, unary_force_functors);
 
     auto start_time = std::chrono::high_resolution_clock::now();
     for (size_t n = 0; n < n_steps; n ++) {