influenceFn.h

//  Copyright (c) 2019 Prashant K. Jha
//  Copyright (c) 2019 Patrick Diehl
//
//  Distributed under the Boost Software License, Version 1.0. (See accompanying
//  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
 
#ifndef MATERIAL_PD_INFLUENCEFN_H
#define MATERIAL_PD_INFLUENCEFN_H
 
#include <cstring>
#include <vector>
#include <cmath>
#include "util/utilIO.h"
 
namespace material {
 
namespace pd {
 
class BaseInfluenceFn {
 
public:
  BaseInfluenceFn() = default;
 
  virtual double getInfFn(const double &r) const = 0;
 
  virtual double getMoment(const size_t &i) const = 0;
 
  virtual std::string printStr(int nt, int lvl) const {
 
    auto tabS = util::io::getTabS(nt);
    std::ostringstream oss;
    oss << tabS << "------- BaseInfluenceFn --------" << std::endl << std::endl;
    oss << tabS << "Provides abstraction for different influence function "
                   "types" << std::endl;
    oss << tabS << std::endl;
 
    return oss.str();
  }
 
  virtual void print(int nt, int lvl) const { std::cout << printStr(nt, lvl); }
 
  virtual void print() const { print(0, 0); }
};
 
class ConstInfluenceFn : public BaseInfluenceFn {
 
public:
  ConstInfluenceFn(const std::vector<double> &params, const size_t &dim): BaseInfluenceFn(), d_a0(0.) {
 
    d_a0 = params.empty() ? double(dim + 1) : params[0];
  }
 
  double getInfFn(const double &r) const override {
    return d_a0;
  }
 
  double getMoment(const size_t &i) const override {
    return d_a0 / double(i + 1);
  }
 
  std::string printStr(int nt, int lvl) const override {
 
    auto tabS = util::io::getTabS(nt);
    std::ostringstream oss;
    oss << tabS << "------- ConstInfluenceFn --------" << std::endl << std::endl;
    oss << tabS << "Constant function with constant = " << d_a0 << std::endl;
    oss << tabS << "First moment = " << getMoment(1)
        << ", second moment = " << getMoment(2)
        << ", third moment = " << getMoment(3) << std::endl;
    oss << tabS << std::endl;
 
    return oss.str();
  }
 
  void print(int nt, int lvl) const override {
    std::cout << printStr(nt, lvl);
  }
 
  void print() const override { print(0, 0); }
 
private:
  double d_a0;
};
 
class LinearInfluenceFn : public BaseInfluenceFn {
 
public:
  LinearInfluenceFn(const std::vector<double> &params, const size_t &dim) : BaseInfluenceFn(), d_a0(0.), d_a1(0.) {
 
    if (params.empty()) {
      // choose a0, a1 = -a0 such that \int_0^1 J(r) r^d dr = 1
      // and J(r) = a0 (1 - r)
      if (dim == 1) {
        d_a0 = 6.;
        d_a1 = -d_a0;
      } else if (dim == 2) {
        d_a0 = 12.;
        d_a1 = -d_a0;
      } else if (dim == 3) {
        d_a0 = 20.;
        d_a1 = -d_a0;
      }
    } else {
      d_a0 = params[0];
      if (params.size() < 2)
        d_a1 = -d_a0;
      else
        d_a1 = params[1];
    }
  }
 
  double getInfFn(const double &r) const override {
    return d_a0 + d_a1 * r;
  }
 
  double getMoment(const size_t &i) const override {
    return (d_a0 / double(i + 1)) + (d_a1 / double(i + 2));
  }
 
  std::string printStr(int nt, int lvl) const override {
 
    auto tabS = util::io::getTabS(nt);
    std::ostringstream oss;
    oss << tabS << "------- LinearInfluenceFn --------" << std::endl << std::endl;
    oss << tabS << "Linear function a0 + a1*r with constants: a0 = "
        << d_a0 << ", a1 = " << d_a1 << std::endl;
    oss << tabS << "First moment = " << getMoment(1)
        << ", second moment = " << getMoment(2)
        << ", third moment = " << getMoment(3) << std::endl;
    oss << tabS << std::endl;
 
    return oss.str();
  }
 
  void print(int nt, int lvl) const override {
    std::cout << printStr(nt, lvl);
  }
 
  void print() const override { print(0, 0); }
 
private:
  double d_a0;
 
  double d_a1;
};
 
class GaussianInfluenceFn : public BaseInfluenceFn {
 
public:
  GaussianInfluenceFn(const std::vector<double> &params, const size_t &dim) : BaseInfluenceFn(), d_alpha(0.), d_beta(0.) {
 
    if (params.empty()) {
      // beta = 0.2 (default value)
      // choose alpha such that \int_0^1 J(r) r^d dr = 1
      d_beta = 0.2;
      if (dim == 1)
        d_alpha = 2. / (d_beta * (1. - std::exp(-1. / d_beta)));
      else if (dim == 2)
        d_alpha = (4.0 / d_beta) * 1.0 /
                  (std::sqrt(M_PI * d_beta) * std::erf(1.0 / std::sqrt(d_beta)) -
                   2.0 * std::exp(-1.0 / d_beta));
      else if (dim == 3)
        d_alpha = (2.0 / d_beta) * 1.0 /
                  (d_beta - (d_beta + 1.) * std::exp(-1.0 / d_beta));
    } else {
      d_alpha = params[0];
      d_beta = params[1];
    }
  }
 
  double getInfFn(const double &r) const override {
    return d_alpha * std::exp(-r * r / d_beta);
  }
 
  double getMoment(const size_t &i) const override {
 
    double sq1 = std::sqrt(d_beta);
    double sq2 = std::sqrt(M_PI);
    // M_i = \int_0^1 alpha exp(-r^2/beta) r^i dr
 
    if (i == 0) {
      // M0 = 0.5 * \alpha (\beta)^(1/2) * (pi)^(1/2) * erf((1/beta)^(1/2))
 
      return 0.5 * d_alpha * sq1 * sq2 * std::erf(1. / sq1);
    } else if (i == 1) {
      // M1 = 0.5 * \alpha \beta (1 - exp(-1/beta))
 
      return 0.5 * d_alpha * d_beta * (1. - std::exp(-1. / d_beta));
    } else if (i == 2) {
      // M2 = 0.5 * \alpha (\beta)^(3/2) * [0.5 * (pi)^(1/2) erf((1/beta)^(1/2)
      // ) - (1/beta)^(1/2) * exp(-1/beta) ]
 
      return 0.5 * d_alpha * d_beta * sq1 *
             (0.5 * sq2 * std::erf(1. / sq1) -
              (1. / sq1) * std::exp(-1. / d_beta));
    } else if (i == 3) {
      // M3 = 0.5 * \alpha (\beta)^(2) * [1 - ((1/beta) + 1) * exp(-1/beta)]
 
      return 0.5 * d_alpha * d_beta * d_beta *
             (1. - (1. + 1. / d_beta) * std::exp(-1. / d_beta));
    }
  }
 
  std::string printStr(int nt, int lvl) const override {
 
    auto tabS = util::io::getTabS(nt);
    std::ostringstream oss;
    oss << tabS << "------- GaussianInfluenceFn --------" << std::endl << std::endl;
    oss << tabS << "Gaussian function a0 * exp(-r*r / a1) with constants: a0 = "
        << d_alpha << ", a1 = " << d_beta << std::endl;
    oss << tabS << "First moment = " << getMoment(1)
        << ", second moment = " << getMoment(2)
        << ", third moment = " << getMoment(3) << std::endl;
    oss << tabS << std::endl;
 
    return oss.str();
  }
 
  void print(int nt, int lvl) const override {
    std::cout << printStr(nt, lvl);
  }
 
  void print() const override { print(0, 0); }
 
private:
  double d_alpha;
 
  double d_beta;
};
 
} // namespace pd
 
} // namespace material
 
#endif // MATERIAL_PD_INFLUENCEFN_H