dirichletBoundaryCondition.hpp

#ifndef MARIAN_DIRICHLETBOUNDARYCONDITION_H
#define MARIAN_DIRICHLETBOUNDARYCONDITION_H


#include <vector>
#include <FDM/tridiagonalOperator.hpp>
#include <FDM/boundaryConditions/boundaryCondition.hpp>

namespace marian {

  template<typename F>
  class DirichletBoundaryCondition: public DCBoundaryCondition<DirichletBoundaryCondition<F> > {
  public:
    DirichletBoundaryCondition() {};

    DirichletBoundaryCondition(BCSide side, F value):
      side_(side), value_(value) {};

    void beforeExplicitStep(TridiagonalOperator&) override;
    void afterExplicitStep(std::vector<double>&, double t) override;

    void beforeImplicitStep(TridiagonalOperator&, std::vector<double>&, double t) override;
    void afterImplicitStep(std::vector<double>&, double t) override;

    std::string info() const override {
      std::string side;
      switch(side_) {
      case BCSide::LOW:  side = " low"; break;
      case BCSide::UPP:  side = " up";  break;
      case BCSide::FREE: side = " free";break;  
      }
      return "DirichletBC" + side;
    }
    virtual ~DirichletBoundaryCondition(){};

  private:
    BCSide side_;    
    F value_;         
  };

  
  template<typename F>
  void DirichletBoundaryCondition<F>::beforeExplicitStep(TridiagonalOperator& op) {
    switch (side_) {
    case BCSide::LOW:
      op.setFirstRow(1.0, 0.0);
      break;
    case BCSide::UPP:
      op.setLastRow(0.0, 1.0);
      break;
    case BCSide::FREE:
      break;
    }
  }

  template<typename F>
  void DirichletBoundaryCondition<F>::afterExplicitStep(std::vector<double>& f, double t) {
    switch (side_) {
    case BCSide::LOW:
      f.front() = value_(t);
      break;
    case BCSide::UPP:
      f.back() = value_(t);
      break;
    case BCSide::FREE:
      break;
    }
  }

  template<typename F>
  void DirichletBoundaryCondition<F>::beforeImplicitStep(TridiagonalOperator& L,
                             std::vector<double>& f,
                             double t) {
    switch (side_) {
    case BCSide::LOW:
      L.setFirstRow(1.0, 0.0);
      f.front() = value_(t);
      break;
    case BCSide::UPP:
      L.setLastRow(0.0, 1.0);
      f.back() = value_(t);
      break;
    case BCSide::FREE:
      break;
    }
  }

  template<typename F>
  void DirichletBoundaryCondition<F>::afterImplicitStep(std::vector<double>&, double) {
  }
  
} // namespace marian


#endif /* MARIAN_DIRICHLETBOUNDARYCONDITION_H */