vtkgridfunction.hh

#ifndef VTK_BASIS_GRID_FUNCTION_HH
#define VTK_BASIS_GRID_FUNCTION_HH

#include <dune/grid/io/file/vtk/function.hh>

namespace Dune {
/** \brief A VTK basis grid function.
 *
 *	This function "evaluates" by evaluating the global basis and
 *interpolating the function values.
 *  \tparam Basis   The global basis.
 *  \tparam CoefficientType The vector type for the coefficients.
*/
template <class Basis, class CoefficientType>
class VTKBasisGridFunction
    : public VTKFunction<typename Basis::GridView::Grid> {

  typedef VTKFunction<typename Basis::GridView::Grid> Base;
  typedef typename Basis::LocalFiniteElement::Traits::LocalBasisType::Traits::
      RangeType RangeType;

public:
  typedef typename Base::Entity Entity;
  typedef typename Base::ctype ctype;
  using Base::dim;

  /** \brief Construct from given global basis, coefficient vector and name.
   *
   *  \param basis    The global basis.
   *  \param v    A corresponding vector of coefficients.
   *  \param s    A name of the function.
   */
  VTKBasisGridFunction(const Basis& basis, const CoefficientType& v,
                       const std::string& s)
      : basis_(basis), coeffs_(v), s(s_) {
    if (v.size() != basis_.size())
      DUNE_THROW(
          IOError,
          "VTKGridFunction: Coefficient vector is not matching the basis");
  }

  /** \brief Get the number of components the function has. */
  virtual int ncomps() const { return CoefficientType::block_type::dimension; }

  /** \brief Locally evaluate a component of the function.
   *
   *  \param comp The component to evaluate.
   *  \param e    The element the local coordinates are taken from.
   *  \param xi   The local coordinates where to evaluate the function.
   */
  virtual double evaluate(int comp, const Entity& e,
                          const Dune::FieldVector<ctype, dim>& xi) const {
    // evaluate the local shapefunctions
    const typename Basis::LocalFiniteElement& localFE =
        basis_.getLocalFiniteElement(e);
    std::vector<RangeType> values(localFE.localBasis().size());
    localFE.localBasis().evaluateFunction(xi, values);

    // get the comp'th component
    double rt = 0;
    for (int i = 0; i < localFE.localBasis().size(); i++)
      rt += values[i] * coeffs_[basis_.index(e, i)][comp];

    return rt;
  }

  /** \brief Stupid function to return the name. */
  virtual std::string name() const { return s_; }

  /** \brief Destructor. */
  virtual ~VTKBasisGridFunction() {}

private:
  const Basis& basis_;
  const CoefficientType& coeffs_;
  std::string s_;
};
}
#endif