gravity.hh

#ifndef DUNE_TECTONIC_GRAVITY_HH
#define DUNE_TECTONIC_GRAVITY_HH

#include <dune/common/function.hh>
#include <dune/common/fvector.hh>

template <int dimension>
class Gravity
    : public Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
                                   Dune::FieldVector<double, dimension>> {
public:
  Gravity(
      Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
                            Dune::FieldVector<double, 1>> const &_densityField,
      Dune::FieldVector<double, dimension> const &_zenith, double _gravity)
      : densityField(_densityField), zenith(_zenith), gravity(_gravity) {}

  void evaluate(Dune::FieldVector<double, dimension> const &x,
                Dune::FieldVector<double, dimension> &y) const {
    y = zenith;

    Dune::FieldVector<double, 1> density;
    densityField.evaluate(x, density);

    y *= -gravity * density;
  }

private:
  Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
                        Dune::FieldVector<double, 1>> const &densityField;
  Dune::FieldVector<double, dimension> const &zenith;
  double const gravity;
};
#endif