bem/kernel.hh

Go to the documentation of this file.

/* Kernel functions for integral equations
 */

#ifndef bemKernel_hh
#define bemKernel_hh

#include "basics/vectorsMatrices.hh"

#include <cmath>
#ifndef M_1_PI
#  define M_1_PI         0.31830988618379067154  /* 1/pi */
#endif

namespace bem {

#define KrnlEPS 1e-10

  // ****************************************************************** Cnst **

  class Cnst {
    friend std::ostream& operator<<(std::ostream& os, const Cnst& c);
    concepts::Real c_;

  public:
    typedef concepts::Real F;

    Cnst(concepts::Real c) : c_(c) {}

    concepts::Real operator()(const concepts::Real3d&, const concepts::Real3d&)
      { return c_; }
  };

  // ****************************************************************** Poly **

  class Poly {
    friend std::ostream& operator<<(std::ostream& os, const Poly& p);
    bool var_;
    uint idx_;
    uint exp_;

  public:
    typedef concepts::Real F;

    Poly(bool var, uint idx, uint e)
      : var_(var), idx_((idx > 2 ? 2 : idx)), exp_(e) {}

    concepts::Real operator()(const concepts::Real3d& x,
            const concepts::Real3d& y) {
      concepts::Real vari = var_ ? y[idx_] : x[idx_];
      concepts::Real krnl = 1.0;
      for(uint i = 0; i < exp_; i++)  krnl *= vari;
      return krnl;
    }
  };

  // ***************************************************************** ExpR2 **

  class ExpR2 {
    friend std::ostream& operator<<(std::ostream& os, const ExpR2& e);
    concepts::MapReal3d gamma_;

  public:
    typedef concepts::Real F;

    ExpR2(concepts::MapReal3d& gamma) : gamma_(gamma) {}

    concepts::Real operator()(const concepts::Real3d& x,
            const concepts::Real3d& y) {
      concepts::Real3d r(x-y);
      concepts::Real   e = r * (gamma_ * r);

      return std::exp(-e);
    }
  };

  // ****************************************************************** ExpR **

  class ExpR {
    friend std::ostream& operator<<(std::ostream& os, const ExpR& e);

    concepts::Real gamma_;

  public:
    typedef concepts::Real F;
    ExpR(concepts::Real gamma) : gamma_(gamma) {}

    concepts::Real operator()(const concepts::Real3d& x,
            const concepts::Real3d& y) {
      concepts::Real r = (x - y).l2();
      return std::exp(-gamma_ * r);
    }
  };

  // **************************************************************** Inv1Rn **

  class Inv1Rn {
    friend std::ostream& operator<<(std::ostream& os, const Inv1Rn& krnl);
    concepts::Real gamma_;
    uint n_;

  public:
    typedef concepts::Real F;

    Inv1Rn(concepts::Real gamma, uint n) : gamma_(gamma), n_(n) {}

    concepts::Real operator()(const concepts::Real3d& x,
            const concepts::Real3d& y) {
      concepts::Real r = (x - y).l2();
      return (r < KrnlEPS) ? 1.0
  : 1.0 / (1.0 + gamma_*std::exp(n_ * std::log(r)));
    }
  };

  // *************************************************************** Laplace **

  class Laplace {
    friend std::ostream& operator<<(std::ostream& os, const Laplace& lpl);
    const concepts::Real sc_;

  public:
    typedef concepts::Real F;

    Laplace() : sc_(0.25 * M_1_PI) {}

    concepts::Real operator()(const concepts::Real3d& x,
            const concepts::Real3d& y) {
      concepts::Real r = (x - y).l2();
      return sc_ / r;
    }
  };

} // namespace bem

#endif // bemKernel_hh