From 9bcccb857e9d26d9da72c0f6915fcb5f68d8a2d0 Mon Sep 17 00:00:00 2001 From: Egor Demidov Date: Sun, 28 Jan 2024 13:59:03 -0500 Subject: [PATCH] Adapted to the changes in libtimestep --- deps/libtimestep | 2 +- .../libgran/granular_system/granular_system.h | 87 ++++++++++++++++--- .../libgran/no_unary_force/no_unary_force.h | 26 ++++++ main.cpp | 24 ++++- test/contact.cpp | 14 ++- test/hamaker.cpp | 22 ++++- test/hamaker_2.cpp | 22 ++++- test/hamaker_3.cpp | 22 ++++- test/sintered.cpp | 17 +++- 9 files changed, 205 insertions(+), 31 deletions(-) create mode 100644 include/libgran/no_unary_force/no_unary_force.h diff --git a/deps/libtimestep b/deps/libtimestep index 22d2bb3..28cbfa6 160000 --- a/deps/libtimestep +++ b/deps/libtimestep @@ -1 +1 @@ -Subproject commit 22d2bb3f92d340f3af9d621c78cf6e55241c9af3 +Subproject commit 28cbfa675c4806031e882b1912ac4421198047a2 diff --git a/include/libgran/granular_system/granular_system.h b/include/libgran/granular_system/granular_system.h index e7974bf..8ab4345 100644 --- a/include/libgran/granular_system/granular_system.h +++ b/include/libgran/granular_system/granular_system.h @@ -14,6 +14,60 @@ std::pair operator + (std::pair const & lhs, std::pair 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_container_t; + + binary_force_functor_container(binary_force_functors_t & ... force_functors) + : binary_force_functors{force_functors...} {} + + std::pair 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::value > 0, "at least one binary force functor must be provided " + "as a template parameter"); + + std::pair accelerations = std::apply([i, j, &x, &v, &theta, &omega, t] (auto & ... e) -> std::pair { + return (e(i, j, x, v, theta, omega, t) + ...); + }, binary_force_functors); + return accelerations; + } + + std::tuple 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_container_t; + + unary_force_functor_container(unary_force_functors_t & ... force_functors) + : unary_force_functors{force_functors...} {} + + std::pair 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::value > 0, "at least one binary force functor must be provided " + "as a template parameter"); + + std::pair accelerations = std::apply([i, &x, &v, &theta, &omega, t] (auto & ... e) -> std::pair { + return (e(i, x, v, theta, omega, t) + ...); + }, unary_force_functors); + return accelerations; + } + + std::tuple 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> { + granular_system> { public: typedef std::vector field_container_t; @@ -42,28 +98,35 @@ class granular_system : public rotational_binary_system & 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> + granular_system> (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 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::value > 0, "at least one force functor must be provided " - "as a template parameter"); - std::pair accelerations = std::apply([i, j, &x, &v, &theta, &omega, t] (auto & ... e) -> std::pair { - 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 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; + binary_force_functor_container_t binary_force_functors; + unary_force_functor_container_t unary_force_functors; }; #endif //LIBGRAN_GRANULAR_SYSTEM_H diff --git a/include/libgran/no_unary_force/no_unary_force.h b/include/libgran/no_unary_force/no_unary_force.h new file mode 100644 index 0000000..faf2ac5 --- /dev/null +++ b/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 +struct no_unary_force_functor { + explicit no_unary_force_functor(field_value_t zero_field) : zero_field(std::move(zero_field)) {} + + std::pair operator () (size_t i [[maybe_unused]], + std::vector const & x [[maybe_unused]], + std::vector const & v [[maybe_unused]], + std::vector const & theta [[maybe_unused]], + std::vector 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 diff --git a/main.cpp b/main.cpp index f9c2a72..d37d49f 100644 --- a/main.cpp +++ b/main.cpp @@ -10,6 +10,7 @@ #include #include +#include #include #include "writer.h" @@ -100,14 +101,29 @@ int main() { hamaker_functor hamaker_model(A, h0, r_part, mass, Eigen::Vector3d::Zero(), 0.0); + binary_force_functor_container, + hamaker_functor> + binary_force_functors{contact_force_model, hamaker_model}; + + no_unary_force_functor no_unary_force(Eigen::Vector3d::Zero()); + + unary_force_functor_container> + 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, - hamaker_functor> 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, + hamaker_functor>, + unary_force_functor_container>> 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 diff --git a/test/contact.cpp b/test/contact.cpp index a0c73f5..5a73072 100644 --- a/test/contact.cpp +++ b/test/contact.cpp @@ -9,6 +9,7 @@ #include #include +#include #include #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> + binary_force_functors{contact_force_model}; + + no_unary_force_functor no_unary_force(Eigen::Vector3d::Zero()); + + unary_force_functor_container> + 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> 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>, + unary_force_functor_container>> 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 ++) { diff --git a/test/hamaker.cpp b/test/hamaker.cpp index 3582e69..8349976 100644 --- a/test/hamaker.cpp +++ b/test/hamaker.cpp @@ -10,6 +10,7 @@ #include #include +#include #include #include "compute_energy.h" @@ -83,14 +84,29 @@ int main() { hamaker_functor hamaker_model(A, h0, r_part, mass, Eigen::Vector3d::Zero(), 0.0); + binary_force_functor_container, + hamaker_functor> + binary_force_functors{contact_force_model, hamaker_model}; + + no_unary_force_functor no_unary_force(Eigen::Vector3d::Zero()); + + unary_force_functor_container> + 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, - hamaker_functor> system(x0, + rotational_step_handler, + binary_force_functor_container, + hamaker_functor>, + unary_force_functor_container>> 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 ++) { diff --git a/test/hamaker_2.cpp b/test/hamaker_2.cpp index 8654dc7..8365905 100644 --- a/test/hamaker_2.cpp +++ b/test/hamaker_2.cpp @@ -10,6 +10,7 @@ #include #include +#include #include #include "../writer.h" @@ -91,14 +92,29 @@ int main() { hamaker_functor hamaker_model(A, h0, r_part, mass, Eigen::Vector3d::Zero(), 0.0); + binary_force_functor_container, + hamaker_functor> + binary_force_functors{contact_force_model, hamaker_model}; + + no_unary_force_functor no_unary_force(Eigen::Vector3d::Zero()); + + unary_force_functor_container> + 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, - hamaker_functor> system(x0, + rotational_step_handler, + binary_force_functor_container, + hamaker_functor>, + unary_force_functor_container>> 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 ++) { diff --git a/test/hamaker_3.cpp b/test/hamaker_3.cpp index 6753d9e..6a63ac5 100644 --- a/test/hamaker_3.cpp +++ b/test/hamaker_3.cpp @@ -10,6 +10,7 @@ #include #include +#include #include #include "../writer.h" @@ -97,14 +98,29 @@ int main() { hamaker_functor hamaker_model(A, h0, r_part, mass, Eigen::Vector3d::Zero(), 0.0); + binary_force_functor_container, + hamaker_functor> + binary_force_functors{contact_force_model, hamaker_model}; + + no_unary_force_functor no_unary_force(Eigen::Vector3d::Zero()); + + unary_force_functor_container> + 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, - hamaker_functor> system(x0, + rotational_step_handler, + binary_force_functor_container, + hamaker_functor>, + unary_force_functor_container>> 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 ++) { diff --git a/test/sintered.cpp b/test/sintered.cpp index 445e99b..31978c0 100644 --- a/test/sintered.cpp +++ b/test/sintered.cpp @@ -11,6 +11,7 @@ #include #include +#include #include #include "mass_distribution.h" @@ -92,6 +93,14 @@ int main() { sinter_functor sinter_model(x0.size(), k, gamma_n, 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, x0.begin(), generate_random_unit_vector, 1.0e-9); + binary_force_functor_container> + binary_force_functors{sinter_model}; + + no_unary_force_functor no_unary_force(Eigen::Vector3d::Zero()); + + unary_force_functor_container> + unary_force_functors(no_unary_force); + // We need to use a custom step handler with the sintering model // Get the custom step handler instance from the sinter model auto step_handler_instance = sinter_model.get_step_handler>(); @@ -100,9 +109,11 @@ int main() { // Using velocity Verlet integrator for rotational systems and a default // step handler for rotational systems granular_system> system(x0, - v0, theta0, omega0, 0.0, Eigen::Vector3d::Zero(), 0.0, step_handler_instance, - sinter_model); + sinter_step_handler_double, binary_force_functor_container>, + unary_force_functor_container>> system(x0, + v0, theta0, omega0, 0.0, Eigen::Vector3d::Zero(), + 0.0, step_handler_instance, + binary_force_functors, unary_force_functors); Eigen::Vector3d center_0 = center_of_mass(system.get_x());