Merge #3265

3265: Bit Mask utility function r=RudolfWeeber a=fweik

Description of changes:
 - tuple interface for Utils::{Array, Vector}
 - Extract bit masking utility function to utils
 - Some small code modernization



Co-authored-by: Florian Weik <[email protected]>

src/core/rotation.cpp

@@ -51,6 +51,7 @@
 #include <utils/math/rotation_matrix.hpp>
 
 #include <cmath>
+#include <utils/mask.hpp>
 
 /** Calculate the derivatives of the quaternion and angular acceleration
  *  for a given particle
@@ -135,19 +136,14 @@ void define_Qdd(Particle const &p, double Qd[4], double Qdd[4], double S[3],
 void propagate_omega_quat_particle(Particle &p) {
 
   // If rotation for the particle is disabled entirely, return early.
-  if (!p.p.rotation)
+  if (p.p.rotation == ROTATION_FIXED)
     return;
 
   Utils::Vector4d Qd{}, Qdd{};
   Utils::Vector3d S{}, Wd{};
 
   // Clear rotational velocity for blocked rotation axes.
-  if (!(p.p.rotation & ROTATION_X))
-    p.m.omega[0] = 0;
-  if (!(p.p.rotation & ROTATION_Y))
-    p.m.omega[1] = 0;
-  if (!(p.p.rotation & ROTATION_Z))
-    p.m.omega[2] = 0;
+  p.m.omega = Utils::mask(p.p.rotation, p.m.omega);
 
   define_Qdd(p, Qd.data(), Qdd.data(), S.data(), Wd.data());
 
@@ -172,25 +168,12 @@ void propagate_omega_quat_particle(Particle &p) {
 }
 
 inline void convert_torque_to_body_frame_apply_fix_and_thermostat(Particle &p) {
-  auto const t = convert_vector_space_to_body(p, p.f.torque);
-  p.f.torque = Utils::Vector3d{};
-
-  if (thermo_switch & THERMO_LANGEVIN) {
-    friction_thermo_langevin_rotation(p);
-
-    p.f.torque += t;
-  } else {
-    p.f.torque = t;
-  }
-
-  if (!(p.p.rotation & ROTATION_X))
-    p.f.torque[0] = 0;
-
-  if (!(p.p.rotation & ROTATION_Y))
-    p.f.torque[1] = 0;
+  auto const torque = (thermo_switch & THERMO_LANGEVIN)
+                          ? friction_thermo_langevin_rotation(p) +
+                                convert_vector_space_to_body(p, p.f.torque)
+                          : convert_vector_space_to_body(p, p.f.torque);
 
-  if (!(p.p.rotation & ROTATION_Z))
-    p.f.torque[2] = 0;
+  p.f.torque = mask(p.p.rotation, torque);
 }
 
 /** convert the torques to the body-fixed frames and propagate angular
@@ -290,31 +273,19 @@ Utils::Vector3d convert_vector_space_to_body(const Particle &p,
  */
 void local_rotate_particle(Particle &p, const Utils::Vector3d &axis_space_frame,
                            const double phi) {
-  // Convert rotation axis to body-fixed frame
-  Utils::Vector3d axis = convert_vector_space_to_body(p, axis_space_frame);
-
   // Rotation turned off entirely?
   if (!p.p.rotation)
     return;
 
-  // Per coordinate fixing
-  if (!(p.p.rotation & ROTATION_X))
-    axis[0] = 0;
-  if (!(p.p.rotation & ROTATION_Y))
-    axis[1] = 0;
-  if (!(p.p.rotation & ROTATION_Z))
-    axis[2] = 0;
-  // Re-normalize rotation axis
-  double l = axis.norm();
-  // Check, if the rotation axis is nonzero
-  if (l < std::numeric_limits<double>::epsilon())
-    return;
-
-  axis /= l;
-
-  double s = sin(phi / 2);
-  Utils::Vector4d q = {cos(phi / 2), s * axis[0], s * axis[1], s * axis[2]};
-  q.normalize();
+  // Convert rotation axis to body-fixed frame
+  auto const axis =
+      mask(p.p.rotation, convert_vector_space_to_body(p, axis_space_frame))
+          .normalize();
+
+  auto const s = std::sin(phi / 2);
+  auto const q =
+      Utils::Vector4d{cos(phi / 2), s * axis[0], s * axis[1], s * axis[2]}
+          .normalize();
 
   // Rotate the particle
   p.r.quat = Utils::multiply_quaternions(p.r.quat, q);

src/core/thermostat.hpp

@@ -261,7 +261,7 @@ inline Utils::Vector3d friction_thermo_langevin(Particle const &p) {
    \xi_i\f$.
     The same friction coefficient \f$\gamma\f$ is used as that for translation.
 */
-inline void friction_thermo_langevin_rotation(Particle &p) {
+inline Utils::Vector3d friction_thermo_langevin_rotation(const Particle &p) {
   extern Thermostat::GammaType langevin_pref2_rotation;
   Thermostat::GammaType langevin_pref_friction_buf, langevin_pref_noise_buf;
 
@@ -297,28 +297,15 @@ inline void friction_thermo_langevin_rotation(Particle &p) {
   }
 #endif /* LANGEVIN_PER_PARTICLE */
 
-  // Rotational degrees of virtual sites are thermostatted,
-  // so no switching here
-
   // Here the thermostats happens
-  Utils::Vector3d noise = v_noise(p.p.identity);
-  for (int j = 0; j < 3; j++) {
+  auto const noise = v_noise(p.p.identity);
 #ifdef PARTICLE_ANISOTROPY
-    if (langevin_pref_noise_buf[j] > 0.0) {
-      p.f.torque[j] = -langevin_pref_friction_buf[j] * p.m.omega[j] +
-                      langevin_pref_noise_buf[j] * noise[j];
-    } else {
-      p.f.torque[j] = -langevin_pref_friction_buf[j] * p.m.omega[j];
-    }
+  return -hadamard_product(langevin_pref_friction_buf, p.m.omega) +
+         hadamard_product(langevin_pref_noise_buf, noise);
 #else
-    if (langevin_pref_noise_buf > 0.0) {
-      p.f.torque[j] = -langevin_pref_friction_buf * p.m.omega[j] +
-                      langevin_pref_noise_buf * noise[j];
-    } else {
-      p.f.torque[j] = -langevin_pref_friction_buf * p.m.omega[j];
-    }
+  return -langevin_pref_friction_buf * p.m.omega +
+         langevin_pref_noise_buf * noise;
 #endif
-  }
 }
 
 #endif // ROTATION

src/utils/include/utils/Array.hpp

@@ -20,6 +20,7 @@
 #define UTILS_ARRAY_HPP
 
 #include "device_qualifier.hpp"
+#include "get.hpp"
 
 #include <boost/serialization/access.hpp>
 
@@ -158,5 +159,19 @@ template <typename T, std::size_t N> struct Array {
     ar &m_storage;
   }
 };
+
+template <std::size_t I, class T, std::size_t N>
+struct tuple_element<I, Array<T, N>> {
+  using type = T;
+};
+
+template <class T, std::size_t N>
+struct tuple_size<Array<T, N>> : std::integral_constant<std::size_t, N> {};
+
+template <std::size_t I, class T, std::size_t N>
+auto get(Array<T, N> const &a) -> std::enable_if_t<(I < N), const T &> {
+  return a[I];
+}
+
 } // namespace Utils
 #endif

src/utils/include/utils/Vector.hpp

@@ -368,5 +368,17 @@ template <typename T, size_t N> struct decay_to_scalar<Vector<T, N>> {
 
 template <typename T> struct decay_to_scalar<Vector<T, 1>> { using type = T; };
 
+template <std::size_t I, class T, std::size_t N>
+struct tuple_element<I, Vector<T, N>> {
+  using type = T;
+};
+
+template <class T, std::size_t N>
+struct tuple_size<Vector<T, N>> : std::integral_constant<std::size_t, N> {};
+
+template <std::size_t I, class T, std::size_t N>
+auto get(Vector<T, N> const &a) -> std::enable_if_t<(I < N), const T &> {
+  return a[I];
+}
 } // namespace Utils
 #endif

src/utils/include/utils/get.hpp

@@ -20,10 +20,16 @@
 #define UTILS_GET_HPP
 
 namespace Utils {
-
 template <std::size_t I, typename T> auto get(const T &v) {
   return std::get<I>(v);
 }
 
+template <class T> struct tuple_size : std::tuple_size<T> {};
+
+template <std::size_t I, class Tuple>
+struct tuple_element : std::tuple_element<I, Tuple> {};
+
+template <std::size_t I, class Tuple>
+using tuple_element_t = typename tuple_element<I, Tuple>::type;
 } // namespace Utils
 #endif

src/utils/include/utils/mask.hpp

@@ -0,0 +1,47 @@
+#ifndef ESPRESSO_MASK_HPP
+#define ESPRESSO_MASK_HPP
+
+#include <utils/Array.hpp>
+#include <utils/get.hpp>
+#include <utils/type_traits.hpp>
+
+#include <utility>
+
+namespace Utils {
+namespace detail {
+template <class T, class Integral, size_t... I>
+auto mask_impl(Integral mask, T t, std::index_sequence<I...>) {
+  return T{((mask & (1u << I)) ? get<I>(t) : tuple_element_t<I, T>{})...};
+}
+} // namespace detail
+
+/**
+ * @brief Pick elements of a tuple-like by a bit mask.
+ *
+ * E.g. every element of the input for which the corresponding
+ * bit is set in the mask is set is copied to the output unmodified,
+ * the elements that are not set are set to zero (default constructed
+ * instance of the type).
+ *
+ * Example:
+ *   mask(0b1011, {1, 2, 3, 4}) => {1, 0, 3, 4}
+ *
+ * @tparam T implements the tuple interface(get, tuple_size, ...)
+ * @tparam Integral An unsigned integral type
+ * @param mask bit mask, if the i-th bit is set, the i-th element
+ *        in @param t is copied to the output, otherwise it is set to zero.
+ * @param t input elements
+ * @return t partially zeroed out according to mask
+ */
+template <class T, class Integral>
+auto mask(Integral mask, T t)
+    -> std::enable_if_t<std::is_unsigned<Integral>::value &&
+                            (size_in_bits<Integral>::value >=
+                             tuple_size<T>::value),
+                        T> {
+  return detail::mask_impl(mask, t,
+                           std::make_index_sequence<tuple_size<T>::value>{});
+}
+} // namespace Utils
+
+#endif // ESPRESSO_MASK_HPP

src/utils/include/utils/type_traits.hpp

@@ -19,6 +19,7 @@
 #ifndef UTILS_TYPE_TRAITS_HPP
 #define UTILS_TYPE_TRAITS_HPP
 
+#include <climits>
 #include <type_traits>
 
 namespace Utils {
@@ -44,6 +45,13 @@ template <class B1, class... Bn>
 struct conjunction<B1, Bn...>
     : std::conditional<bool(B1::value), conjunction<Bn...>, B1>::type {};
 
+/**
+ * @brief Size of a type in bits.
+ *
+ * On posix plattforms this is 8 * sizeof(T).
+ */
+template <class T>
+struct size_in_bits : std::integral_constant<size_t, CHAR_BIT * sizeof(T)> {};
 } // namespace Utils
 
 #endif

src/utils/tests/Array_test.cpp

@@ -103,3 +103,13 @@ BOOST_AUTO_TEST_CASE(zero_size) {
   Array<int, 0> const a{};
   BOOST_CHECK(a.empty());
 }
+
+BOOST_AUTO_TEST_CASE(tuple_protocol) {
+  using A = Utils::Array<int, 4>;
+
+  static_assert(std::is_same<Utils::tuple_element_t<0, A>, int>::value, "");
+  static_assert(std::is_same<Utils::tuple_element_t<1, A>, int>::value, "");
+  static_assert(A{}.size() == Utils::tuple_size<A>::value, "");
+
+  BOOST_CHECK_EQUAL(Utils::get<1>(A{1, 2, 3, 4}), 2);
+}

src/utils/tests/CMakeLists.txt

@@ -37,6 +37,8 @@ unit_test(NAME sampling_test SRC sampling_test.cpp DEPENDS utils)
 unit_test(NAME coordinate_transformation_test SRC coordinate_transformation.cpp DEPENDS utils)
 unit_test(NAME rotation_matrix_test SRC rotation_matrix_test.cpp DEPENDS utils)
 unit_test(NAME quaternion_test SRC quaternion_test.cpp DEPENDS utils)
+unit_test(NAME mask_test SRC mask_test.cpp DEPENDS utils)
+unit_test(NAME type_traits_test SRC type_traits_test.cpp DEPENDS utils)
 
 unit_test(NAME gather_buffer_test SRC gather_buffer_test.cpp DEPENDS utils Boost::mpi MPI::MPI_CXX NUM_PROC 4)
 unit_test(NAME scatter_buffer_test SRC scatter_buffer_test.cpp DEPENDS utils Boost::mpi MPI::MPI_CXX NUM_PROC 4)

src/utils/tests/mask_test.cpp

@@ -0,0 +1,47 @@
+/*
+ * Copyright (C) 2019 The ESPResSo project
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#define BOOST_TEST_MODULE mask test
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+#include <utils/mask.hpp>
+
+#include <array>
+#include <cstdint>
+#include <tuple>
+
+BOOST_AUTO_TEST_CASE(mask_) {
+  using Utils::get;
+  using Utils::tuple_element_t;
+
+  const uint8_t mask = 1u | 4u;
+
+  auto const a = std::make_tuple(std::string("a"), 3, 4.5);
+  using input_type = decltype(a);
+
+  auto const result = Utils::mask(mask, a);
+  using result_type = decltype(result);
+
+  static_assert(std::is_same<input_type, result_type>::value, "");
+
+  BOOST_CHECK_EQUAL(get<0>(result), get<0>(a));
+  BOOST_CHECK_EQUAL(get<1>(result), 0);
+  BOOST_CHECK_EQUAL(get<2>(result), get<2>(a));
+}

src/utils/tests/type_traits_test.cpp

@@ -0,0 +1,36 @@
+/*
+ * Copyright (C) 2018-2019 The ESPResSo project
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/** \file
+ * Unit test for Utils Type Traits.
+ *
+ */
+
+#define BOOST_TEST_MODULE type traits tests
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+#include <utils/type_traits.hpp>
+
+#include <climits>
+
+BOOST_AUTO_TEST_CASE(size_in_bits) {
+  static_assert(CHAR_BIT == Utils::size_in_bits<char>::value, "");
+  static_assert(CHAR_BIT * sizeof(int) == Utils::size_in_bits<int>::value, "");
+}