/******************************************************************************
*                                                                             *
*  Copyright 2019 Max Planck Institute for Dynamics and Self-Organization     *
*                                                                             *
*  This file is part of TurTLE.                                               *
*                                                                             *
*  TurTLE 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.                                     *
*                                                                             *
*  TurTLE 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 TurTLE.  If not, see <http://www.gnu.org/licenses/>             *
*                                                                             *
* Contact: Cristian.Lalescu@ds.mpg.de                                         *
*                                                                             *
******************************************************************************/



#ifndef P2P_COMPUTER_HPP
#define P2P_COMPUTER_HPP

#include <cstring>
#include <cassert>

template <class real_number, class partsize_t>
class p2p_computer{

    bool isActive;

    /** \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.
     */
    static double dumb_distance_weight(
            const double dist_pow2,
            const double cutoff){
        // this function should only be called for interacting particles,
        // and particles interact if they are closer than cutoff.
        assert(dist_pow2 < cutoff*cutoff);
        return 1.0;
    }


public:
    p2p_computer() : isActive(true){}

    template <int size_particle_rhs>
    void init_result_array(real_number rhs[], const partsize_t nbParticles) const{
        memset(rhs, 0, sizeof(real_number)*nbParticles*size_particle_rhs);
    }

    template <int size_particle_rhs>
    void reduce_particles_rhs(real_number rhs_dst[], const real_number rhs_src[], const partsize_t nbParticles) const{
        static_assert(size_particle_rhs == 6, "This kernel works only with 6 values per particle's rhs");
        for(int idx_part = 0 ; idx_part < nbParticles ; ++idx_part){
            // We merge only the values modified by the current kernel (3-5)
            for(int idx_rhs = 3 ; idx_rhs < size_particle_rhs ; ++idx_rhs){
                rhs_dst[idx_part*size_particle_rhs+idx_rhs] += rhs_src[idx_part*size_particle_rhs+idx_rhs];
            }
        }
    }

    template <int size_particle_positions, int size_particle_rhs>
    void compute_interaction(const partsize_t /*idx_part1*/,
                             const real_number pos_part1[], real_number rhs_part1[],
                             const partsize_t /*idx_part2*/,
                             const real_number pos_part2[], real_number rhs_part2[],
                             const real_number dist_pow2, const real_number cutoff,
                             const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/) const{
        static_assert(size_particle_positions == 6, "This kernel works only with 6 values for one particle's position+orientation");
        static_assert(size_particle_rhs == 6, "This kernel works only with 6 values per particle's rhs");

        // 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.
        const double ww = dumb_distance_weight(dist_pow2, cutoff);
        ///
        /// term in equation is:
        ///
        /// \f[
        ///     (4 / \tau) \sum_j W_\ell ( | x^i - x^j | ) (p^i \cdot p^j)p^j
        /// \f]
        ///
        const double dot_product = (pos_part1[3+IDXC_X]*pos_part2[3+IDXC_X] +
                              pos_part1[3+IDXC_Y]*pos_part2[3+IDXC_Y] +
                              pos_part1[3+IDXC_Z]*pos_part2[3+IDXC_Z]);
        rhs_part1[3+IDXC_X] += pos_part2[3+IDXC_X] * 4 * ww * dot_product;
        rhs_part1[3+IDXC_Y] += pos_part2[3+IDXC_Y] * 4 * ww * dot_product;
        rhs_part1[3+IDXC_Z] += pos_part2[3+IDXC_Z] * 4 * ww * dot_product;
        rhs_part2[3+IDXC_X] += pos_part1[3+IDXC_X] * 4 * ww * dot_product;
        rhs_part2[3+IDXC_Y] += pos_part1[3+IDXC_Y] * 4 * ww * dot_product;
        rhs_part2[3+IDXC_Z] += pos_part1[3+IDXC_Z] * 4 * ww * dot_product;
    }


    void merge(const p2p_computer&){}

    bool isEnable() const {
        return isActive;
    }

    void setEnable(const bool inIsActive) {
        isActive = inIsActive;
    }
};

#endif