add preliminary particle orientation sync

bfps/cpp/particles/p2p_computer.hpp

@@ -3,6 +3,22 @@
 
 #include <cstring>
 
+/** \brief A simple distance weighting function.
+ *
+ *  This function returns 1 if a distance is smaller than a cut-off length,
+ *  i.e. particle 1 interacts with particle 2 if particle 2 is inside a
+ *  sphere of radius `cutoff' centered on particle 1.
+ */
+double dumb_distance_weight(
+        double distance,
+        double cutoff)
+{
+    // this function should only be called for interacting particles,
+    // and particles interact if they are closer than cutoff.
+    assert(distance < cutoff);
+    return 1.0;
+}
+
 template <class real_number, class partsize_t>
 class p2p_computer{
 public:
@@ -27,10 +43,41 @@ public:
                              const real_number pos_part2[], const real_number data_part2[], real_number rhs_part2[],
                              const real_number dist_pow2,
                              const real_number xshift_coef, const real_number yshift_coef, const real_number zshift_coef) const{
-        // TODO put the kernel here
         static_assert(size_particle_positions == 3, "This kernel works only with 3 values for one position");
+        // TODO:
+        // Should I put 0 in the rhs corresponding to the orientations?
+        // Or should I just add the interaction term?
+        // In other words: is this method called after the vorticity and strain terms have been computed?
+        // The following two lines set the rhs to 0:
+        //std::fill_n(rhs_part1+3, 3, 0);
+        //std::fill_n(rhs_part2+3, 3, 0);
+        real_number distance = sqrt(dist_pow2);
+        double max_distance = 1.0;
+        double tau = 1.0;
+        if (distance >= max_distance)
+            return;
+        // TODO: a reasonable way of choosing between different distance_weight functions should be thought of.
+        // We need to ask Michael about how flexible this distance_weight needs to be.
+        double ww = dumb_distance_weight(distance, max_distance);
+        ///
+        /// term in equation is:
+        ///
+        /// \f[
+        ///     (4 / \tau) \sum_j W_\ell ( | x^i - x^j | ) (p^i \cdot p^j)p^j
+        /// \f]
+        ///
+        double dot_product = (data_part1[0]*data_part2[0] +
+                              data_part1[1]*data_part2[1] +
+                              data_part1[2]*data_part2[2])
+        rhs_part1[3] += data_part2[0] * 4 * ww * dot_product;
+        rhs_part1[4] += data_part2[1] * 4 * ww * dot_product;
+        rhs_part1[5] += data_part2[2] * 4 * ww * dot_product;
+        rhs_part2[3] += data_part1[0] * 4 * ww * dot_product;
+        rhs_part2[4] += data_part1[1] * 4 * ww * dot_product;
+        rhs_part2[5] += data_part1[2] * 4 * ww * dot_product;
     }
 
+
     constexpr static bool isEmpty() {
         return false;
     }