source: src/units/particle/bspline.hpp@ cd0fed

/**
 * @file   bspline.hpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Mon Nov 21 13:27:22 2011
 *
 * @brief  B-Splines for molecular dynamics.
 *
 */

#ifndef BSPLINE_HPP_
#define BSPLINE_HPP_

#include "base/helper.hpp"
#include "base/index.hpp"
#include "base/polynomial.hpp"
#include "base/vector.hpp"
#include "grid/grid.hpp"
#include "units/particle/particle.hpp"

namespace VMG
{

namespace Particle
{

class BSpline
{
public:
  BSpline(const vmg_float& width);

  vmg_float EvaluateSpline(const vmg_float& val) const
  {
    for (unsigned int i=0; i<intervals.size(); ++i)
      if (val < intervals[i])
        return spline_nom[i](val) / spline_denom[i](val);
    return 0.0;
  }

  void SetSpline(Grid& grid, const Particle& p, const int& near_field_cells) const
  {
    assert(p.Pos().X() >= grid.Extent().Begin().X() && p.Pos().X() < grid.Extent().End().X());
    assert(p.Pos().Y() >= grid.Extent().Begin().Y() && p.Pos().Y() < grid.Extent().End().Y());
    assert(p.Pos().Z() >= grid.Extent().Begin().Z() && p.Pos().Z() < grid.Extent().End().Z());

    std::vector<vmg_float> vals(Helper::intpow(2*near_field_cells+1,3));
    vmg_float temp_val;
    vmg_float int_val = 0.0;
    int c = 0;

    Index i;
    Vector dir;

    const Index index_global = static_cast<Index>((p.Pos() - grid.Extent().Begin()) / grid.Extent().MeshWidth());
    assert(index_global.X() >= grid.Global().LocalBegin().X() && index_global.X() < grid.Global().LocalEnd().X());
    assert(index_global.Y() >= grid.Global().LocalBegin().Y() && index_global.Y() < grid.Global().LocalEnd().Y());
    assert(index_global.Z() >= grid.Global().LocalBegin().Z() && index_global.Z() < grid.Global().LocalEnd().Z());

    const Index index_local = index_global - grid.Global().LocalBegin() + grid.Local().Begin();
    assert(index_local.X() >= grid.Local().Begin().X() && index_local.X() < grid.Local().End().X());
    assert(index_local.Y() >= grid.Local().Begin().Y() && index_local.Y() < grid.Local().End().Y());
    assert(index_local.Z() >= grid.Local().Begin().Z() && index_local.Z() < grid.Local().End().Z());

    const Vector pos_beg = p.Pos() - grid.Extent().Begin() - grid.Extent().MeshWidth() * (index_global - near_field_cells);
    const Vector& h = grid.Extent().MeshWidth();

    // Iterate over all grid points which lie in the support of the interpolating B-Spline
    dir[0] = pos_beg[0];
    for (i.X()=-1*near_field_cells; i.X()<=near_field_cells; ++i.X()) {
      dir[1] = pos_beg[1];
      for (i.Y()=-1*near_field_cells; i.Y()<=near_field_cells; ++i.Y()) {
        dir[2] = pos_beg[2];
        for (i.Z()=-1*near_field_cells; i.Z()<=near_field_cells; ++i.Z()) {

          // Compute distance from grid point to particle
          temp_val = EvaluateSpline(dir.Length());
          vals[c++] = temp_val * p.Charge();
          int_val += temp_val;

          dir[2] -= h[2];
        }
        dir[1] -= h[1];
      }
      dir[0] -= h[0];
    }

    // Reciprocal value of the numerically integrated spline
    int_val = 1.0 / (int_val * grid.Extent().MeshWidth().Product());

    c = 0;
    for (i.X()=-1*near_field_cells; i.X()<=near_field_cells; ++i.X())
      for (i.Y()=-1*near_field_cells; i.Y()<=near_field_cells; ++i.Y())
        for (i.Z()=-1*near_field_cells; i.Z()<=near_field_cells; ++i.Z())
          grid(index_local+i) += vals[c++] * int_val;

  }

  vmg_float EvaluatePotential(const vmg_float& val) const
  {
    for (unsigned int i=0; i<intervals.size(); ++i)
      if (val < intervals[i])
        return potential_nom[i](val) / potential_denom[i](val);
    return potential_nom.back()(val) / potential_denom.back()(val);
  }

  vmg_float EvaluateField(const vmg_float& val) const
  {
    for (unsigned int i=0; i<intervals.size(); ++i)
      if (val < intervals[i])
        return field_nom[i](val) / field_denom[i](val);
    return 0.0;
  }

  const vmg_float& GetAntiDerivativeAtZero() const
  {
    return antid;
  }

private:
  std::vector<Polynomial> spline_nom, spline_denom;
  std::vector<Polynomial> potential_nom, potential_denom;
  std::vector<Polynomial> field_nom, field_denom;
  vmg_float antid;
  std::vector<vmg_float> intervals;

  vmg_float R;
};

}

}

#endif /* BSPLINE_HPP_ */