SolidHarmonicsExpansion.h

/*
 * SolidHarmonicsExpansion.h
 *
 *  Created on: Nov 21, 2014
 *      Author: tchipevn, rajats
 */
 
#ifndef SOLIDHARMONICSEXPANSION_H_
#define SOLIDHARMONICSEXPANSION_H_
 
#include "bhfmm/utils/Vector3.h"
#include "bhfmm/expansions/SolidHarmonicsStorage.h"
#include "bhfmm/utils/WignerMatrix.h"
#ifdef FMM_FFT
#include "bhfmm/fft/FFTAccelerableExpansion.h"
#endif  /* FMM_FFT */
#include <cstdlib>
#include <cmath>
 
namespace bhfmm {
class SolidHarmonicsExpansion;
 
void swap(SolidHarmonicsExpansion & s1, SolidHarmonicsExpansion & s2);
 
double potentialML(const SolidHarmonicsExpansion & ME, const SolidHarmonicsExpansion & LE);
 
Vector3<double> forceGradLAndM(const SolidHarmonicsExpansion & LE, const SolidHarmonicsExpansion & ME);
 
Vector3<double> forceLAndGradM(const SolidHarmonicsExpansion & LE, const SolidHarmonicsExpansion & ME);
 
class SolidHarmonicsExpansion
#ifdef FMM_FFT
        : public FFTAccelerableExpansion
#endif  /* FMM_FFT */
{
public:
    SolidHarmonicsExpansion(int order, bool initializeToZero = true);
 
    SolidHarmonicsExpansion(const SolidHarmonicsExpansion & rhs);
 
    virtual ~SolidHarmonicsExpansion();
 
    friend void swap(SolidHarmonicsExpansion & s1, SolidHarmonicsExpansion & s2);
 
    // OPERATORS //
    SolidHarmonicsExpansion & operator =(SolidHarmonicsExpansion RHS);
 
    SolidHarmonicsExpansion & operator+=(const SolidHarmonicsExpansion & RHS);
 
    SolidHarmonicsExpansion & operator*=(double scalar);
 
    // METHODS //
    void clear();
 
    void evaluateLOfR(Vector3<double> r);
 
    void evaluateMOfR(Vector3<double> r);
 
    void convoluteLL(const SolidHarmonicsExpansion & LE1, const SolidHarmonicsExpansion & LE2);
 
    void convoluteLL_Z(const SolidHarmonicsExpansion & LE1, const SolidHarmonicsExpansion & LE2);
 
    void convoluteLM(const SolidHarmonicsExpansion & LE, const SolidHarmonicsExpansion & ME);
 
    void convoluteLM_Z(const SolidHarmonicsExpansion & LE, const SolidHarmonicsExpansion & ME);
 
    void convoluteL_ZM(const SolidHarmonicsExpansion & LE, const SolidHarmonicsExpansion & ME);
 
    void scaleL(double factor);
 
    void scaleM(double factor);
 
    void rotatePhi(const double* CosSinPhi, int negate);
 
    void convoluteWL(const WignerMatrix & W, const SolidHarmonicsExpansion & LE);
 
    void convoluteWM(const WignerMatrix & W, const SolidHarmonicsExpansion & ME);
 
    void clearMonopole();
 
    void setAtMinusR();
 
    friend double potentialML(const SolidHarmonicsExpansion & ME, const SolidHarmonicsExpansion & LE);
 
    friend Vector3<double> forceGradLAndM(const SolidHarmonicsExpansion & LE, const SolidHarmonicsExpansion & ME);
 
    friend Vector3<double> forceLAndGradM(const SolidHarmonicsExpansion & LE, const SolidHarmonicsExpansion & ME);
 
    int getOrder() const {
        return _order;
    }
 
    double & getC(int l, int m) {
        return acc_C(l,m);
    }
  //for FFTAccelerableExpansion
  double & get_C(unsigned l, unsigned m) {
        return acc_C(l,m);
    }
 
    double & getS(int l, int m) {
        return acc_S(l,m);
    }
  //for FFTAccelerableExpansion
  double & get_S(unsigned l, unsigned m) {
        return acc_S(l,m);
    }
 
    int getNumEntries() const {
        return _c.getTotalNumValues() + _s.getTotalNumValues();
    }
 
    void print() const;
 
    void writeValuesToMPIBuffer(double * buf, int& position) const {
        const int end = _c.getTotalNumValues();
        for (int i = 0; i < end; ++i) {
            buf[position++] = acc_c_C_seq(i);
        }
        for (int i = 0; i < end; ++i) {
            buf[position++] = acc_c_S_seq(i);
        }
    }
 
    void readValuesFromMPIBuffer(double * buf, int& position) {
        const int end = _c.getTotalNumValues();
        for (int i = 0; i < end; ++i) {
            acc_C_seq(i) = buf[position++];
        }
        for (int i = 0; i < end; ++i) {
            acc_S_seq(i) = buf[position++];
        }
    }
    void addValuesFromMPIBuffer(double * buf, int& position) {
        const int end = _c.getTotalNumValues();
        for (int i = 0; i < end; ++i) {
            acc_C_seq(i) += buf[position++];
        }
        for (int i = 0; i < end; ++i) {
            acc_S_seq(i) += buf[position++];
        }
    }
 
private:
    //private accessors to C and S terms
 
    double & acc_C(int l, int m) {
        return _c.getValue(l, m);
    }
 
    double & acc_S(int l, int m) {
        return _s.getValue(l, m);
    }
 
    double & acc_C_seq(int r) {
        return _c.getValueSequential(r);
    }
 
    double & acc_S_seq(int r) {
        return _s.getValueSequential(r);
    }
 
    double acc_c_C(int l, int m) const {
        return _c.getValueConst(l, m);
    }
 
    double acc_c_S(int l, int m) const {
        return _s.getValueConst(l, m);
    }
 
    double acc_c_C_seq(int r) const {
        return _c.getValueConstSequential(r);
    }
 
    double acc_c_S_seq(int r) const {
        return _s.getValueConstSequential(r);
    }
 
    void signed_acc_const_CS(int l, int m, double & c, double & s) const {
        const int ind = _c.index(l, abs(m));
        c = _c.getValueConstSequential(ind);
        s = _s.getValueConstSequential(ind);
        if (m >= 0) {
            return;
        } else {
            const int minus_one_pow_m = (m & 1) ? -1.0 : 1.0;
            c *= minus_one_pow_m;
            s *= -minus_one_pow_m;
        }
    }
 
protected:
    int _order;
 
private:
    SolidHarmonicsStorage _c;
 
    SolidHarmonicsStorage _s; // S-terms: imaginary or (-imaginary) part, see 13.4.8
 
};
 
SolidHarmonicsExpansion operator+(SolidHarmonicsExpansion LHS, const SolidHarmonicsExpansion & RHS);
 
SolidHarmonicsExpansion operator*(double scalar, SolidHarmonicsExpansion RHS);
 
SolidHarmonicsExpansion evaluateLOfR(int order, Vector3<double> r);
 
SolidHarmonicsExpansion evaluateMOfR(int order, Vector3<double> r);
 
SolidHarmonicsExpansion convoluteLL(const SolidHarmonicsExpansion & LE1, const SolidHarmonicsExpansion & LE2);
 
SolidHarmonicsExpansion convoluteLL_Z(const SolidHarmonicsExpansion & LE1, const SolidHarmonicsExpansion & LE2);
 
SolidHarmonicsExpansion convoluteLM(const SolidHarmonicsExpansion & LE, const SolidHarmonicsExpansion & ME);
 
SolidHarmonicsExpansion convoluteLM_Z(const SolidHarmonicsExpansion & LE, const SolidHarmonicsExpansion & ME);
 
SolidHarmonicsExpansion convoluteL_ZM(const SolidHarmonicsExpansion & LE, const SolidHarmonicsExpansion & ME);
 
SolidHarmonicsExpansion rotatePhi(const SolidHarmonicsExpansion & E, const double* CosSinPhi, int negate);
 
SolidHarmonicsExpansion rotateThetaL(const SolidHarmonicsExpansion & LE, const WignerMatrix& W);
 
SolidHarmonicsExpansion rotateThetaM(const SolidHarmonicsExpansion & ME, const WignerMatrix& W);
 
SolidHarmonicsExpansion setAtMinusR(SolidHarmonicsExpansion E);
 
} /* namespace bhfmm */
 
#endif /* SOLIDHARMONICSEXPANSION_H_ */