#ifndef SRC_ONE_BODY_PROBLEM_DATA_BC_HH
#define SRC_ONE_BODY_PROBLEM_DATA_BC_HH

#include <cassert>

#include <dune/common/function.hh>

class VelocityDirichletCondition
    : public Dune::VirtualFunction<double, double> {

public:
  VelocityDirichletCondition(const double finalVelocity = 5e-5, const double threshold = 0.1) :
      finalVelocity_(finalVelocity),
      threshold_(threshold)
  {}

  void evaluate(double const &relativeTime, double &y) const {
    // Assumed to vanish at time zero
    
    //std::cout << "VelocityDirichletCondition::evaluate()" << std::endl;
    
    /*if (relativeTime <= 0.1)
        std::cout << "- loading phase" << std::endl;
    else
        std::cout << "- final velocity reached" << std::endl;*/
    
    y = (relativeTime <= threshold_)
            ? finalVelocity_ * (1.0 - std::cos(relativeTime * M_PI / threshold_)) / 2.0
            : finalVelocity_;

    //y = relativeTime * finalVelocity;

  }

private:
  const double finalVelocity_;
  const double threshold_;
};


class VelocityDirichletConditionLinearLoading
    : public Dune::VirtualFunction<double, double> {

public:
  VelocityDirichletConditionLinearLoading(const double finalVelocity = 5e-5, const double threshold = 0.1) :
      finalVelocity_(finalVelocity),
      threshold_(threshold)
  {}

  void evaluate(double const &relativeTime, double &y) const {
    // Assumed to vanish at time zero

    //std::cout << "VelocityDirichletCondition::evaluate()" << std::endl;

    /*if (relativeTime <= 0.1)
        std::cout << "- loading phase" << std::endl;
    else
        std::cout << "- final velocity reached" << std::endl;*/

    y = (relativeTime <= threshold_)
            ? finalVelocity_ * relativeTime / threshold_
            : finalVelocity_;

    //y = relativeTime * finalVelocity;

  }

private:
  const double finalVelocity_;
  const double threshold_;
};

class VelocityStepDirichletCondition
    : public Dune::VirtualFunction<double, double> {

public:
  VelocityStepDirichletCondition(const double initialVelocity = 5e-5, const double finalVelocity = 1e-4, const double loadingThreshold = 0.1, const double stepTime = 0.5) :
      initialVelocity_(initialVelocity),
      finalVelocity_(finalVelocity),
      loadingThreshold_(loadingThreshold),
      stepTime_(stepTime)
  {
      assert(loadingThreshold_ < stepTime_);
  }

  void evaluate(double const &relativeTime, double &y) const {
    // Assumed to vanish at time zero

    y = finalVelocity_;

    if (relativeTime <= loadingThreshold_) {
        y = initialVelocity_ * (1.0 - std::cos(relativeTime * M_PI / loadingThreshold_)) / 2.0;
    } else if (relativeTime < stepTime_) {
        y = initialVelocity_;
    }
  }

private:
  const double initialVelocity_;
  const double finalVelocity_;
  const double loadingThreshold_;
  const double stepTime_;
};

class NeumannCondition : public Dune::VirtualFunction<double, double> {
  public:
    NeumannCondition(double c = 0.0) : c_(c) {}
    void evaluate(double const &relativeTime, double &y) const { y = -c_; }

  private:
    double c_;
};
#endif