forward_euler_problem.h

Go to the documentation of this file.

#ifndef ForwardEulerProblem_H_
#define ForwardEulerProblem_H_

#include "initialproblem.h" 
#include "primal_ts_base.h"

namespace DOpE
{
  template<typename OPTPROBLEM, typename SPARSITYPATTERN, typename VECTOR,
      int dopedim, int dealdim,
      template <int, int> class FE = dealii::FESystem,
      template <int, int> class DH = dealii::DoFHandler>
    class ForwardEulerProblem : public PrimalTSBase<OPTPROBLEM,
    SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>
    {
      public:
        ForwardEulerProblem(OPTPROBLEM& OP) :
            PrimalTSBase<OPTPROBLEM, SPARSITYPATTERN, VECTOR, dopedim, dealdim,
                FE, DH>(OP)
        {
          initial_problem_ = NULL;
        }
        ~ForwardEulerProblem()
        {
          if (initial_problem_ != NULL)
            delete initial_problem_;
        }

      /******************************************************/

      std::string
      GetName()
      {
        return "forward Euler";
      }
      
      /******************************************************/

      InitialProblem<ForwardEulerProblem<OPTPROBLEM, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>, VECTOR, dealdim>&
      GetInitialProblem()
      {
        if (initial_problem_ == NULL)
        {
          initial_problem_ = new InitialProblem<ForwardEulerProblem<OPTPROBLEM, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>, VECTOR, dealdim>
          (*this);
        }
        return *initial_problem_;
      }

      /******************************************************/
      ForwardEulerProblem<OPTPROBLEM, SPARSITYPATTERN, VECTOR, dopedim,
          dealdim, FE, DH>&
      GetBaseProblem()
      {
        return *this;
      }
      /******************************************************/
      
      template<typename EDC>
        void
        ElementEquation(const EDC& edc,
                     dealii::Vector<double> &local_vector, double scale, double /*scale_ico*/)
        {
          if (this->GetPart() == "New")
            {
              dealii::Vector<double> tmp(local_vector);
              tmp = 0.0;
              this->GetProblem().ElementEquation(edc, tmp, 0., scale);
              local_vector += tmp;

              tmp = 0.0;
              this->GetProblem().ElementTimeEquation(edc, tmp,
                  scale );
              local_vector += tmp;

              this->GetProblem().ElementTimeEquationExplicit(edc, local_vector,
                  scale );

            }
          else if (this->GetPart() == "Old")
            {
              dealii::Vector<double> tmp(local_vector);
              tmp = 0.0;
              this->GetProblem().ElementEquation(edc, tmp, scale, 0.);
              local_vector += tmp;

              this->GetProblem().ElementTimeEquation(
                  edc,
                  local_vector,
                  (-1) * scale);
            }
          else
            {
              abort();
            }
        }

      /******************************************************/

     template<typename EDC>
        void
        ElementRhs(const EDC& edc,
            dealii::Vector<double> &local_vector, double scale)
        {
          if (this->GetPart() == "New")
            {

            }
          else if (this->GetPart() == "Old")
            {
              this->GetProblem().ElementRhs(edc, local_vector, scale);
            }
          else
            {
              abort();
            }
        }


        void
        PointRhs(
            const std::map<std::string, const dealii::Vector<double>*> &param_values,
            const std::map<std::string, const VECTOR*> &domain_values,
            VECTOR& rhs_vector, double scale)
        {
          if (this->GetPart() == "New")
          {
          }
          else if (this->GetPart() == "Old")
          {
            this->GetProblem().PointRhs(param_values, domain_values, rhs_vector, scale);
          }
          else
          {
            abort();
          }
        }

      /******************************************************/

      template<typename EDC>
        void
        ElementMatrix(const EDC& edc,
                   dealii::FullMatrix<double> &local_matrix)
        {
          assert(this->GetPart() == "New");
          dealii::FullMatrix<double> m(local_matrix);

          this->GetProblem().ElementMatrix(edc, local_matrix, 0., 1.);

          m = 0.;
          this->GetProblem().ElementTimeMatrix(edc, m);
          local_matrix.add(
              1.0, m);

          m = 0.;
          this->GetProblem().ElementTimeMatrixExplicit(edc, m);
          local_matrix.add(
              1.0 , m);

        }

      /******************************************************/

     template<typename FDC>
        void
        FaceEquation(const FDC& fdc,
                     dealii::Vector<double> &local_vector, double scale, double /*scale_ico*/)
        {
          if (this->GetPart() == "New")
            {
              this->GetProblem().FaceEquation(fdc, local_vector, 0., scale);
            }
          else if (this->GetPart() == "Old")
            {
              this->GetProblem().FaceEquation(fdc, local_vector, scale,0.);
            }
          else
            {
              abort();
            }

        }
      /******************************************************/

      template<typename FDC>
        void
        InterfaceEquation(const FDC& fdc,
                          dealii::Vector<double> &local_vector, double scale, double /*scale_ico*/)
        {
          if (this->GetPart() == "New")
            {
              this->GetProblem().InterfaceEquation(fdc, local_vector, 0., scale);
            }
          else if (this->GetPart() == "Old")
            {
              this->GetProblem().InterfaceEquation(fdc, local_vector, scale,0.);
            }
          else
            {
              abort();
            }

        }

      /******************************************************/

      template<typename FDC>
        void
        FaceRhs(const FDC& fdc,
            dealii::Vector<double> &local_vector, double scale = 1.)
        {
          this->GetProblem().FaceRhs(fdc, local_vector, scale);
        }

      /******************************************************/

      template<typename FDC>
        void
        FaceMatrix(const FDC& fdc,
                   dealii::FullMatrix<double> &local_matrix)
        {
          assert(this->GetPart() == "New");
          this->GetProblem().FaceMatrix(fdc, local_matrix, 0., 1.);
          
        }

      /******************************************************/

      template<typename FDC>
        void
        InterfaceMatrix(const FDC& fdc,
            dealii::FullMatrix<double> &local_matrix)
        {
          assert(this->GetPart() == "New");
          this->GetProblem().InterfaceMatrix(fdc, local_matrix, 0., 1.);
          
        }

      /******************************************************/

      template<typename FDC>
        void
        BoundaryEquation(const FDC& fdc,
                         dealii::Vector<double> &local_vector, double scale, double /*scale_ico*/)
        {
          if (this->GetPart() == "New")
            {
              this->GetProblem().BoundaryEquation(fdc, local_vector, 0., scale);
            }
          else if (this->GetPart() == "Old")
            {
              this->GetProblem().BoundaryEquation(fdc, local_vector, scale,0.);
            }
          else
            {
              abort();
            }

        }

      /******************************************************/

      template<typename FDC>
        void
        BoundaryRhs(const FDC& fdc,
            dealii::Vector<double> &local_vector, double scale)
        {
          this->GetProblem().BoundaryRhs(fdc, local_vector, scale);
        }

      /******************************************************/

      template<typename FDC>
        void
        BoundaryMatrix(const FDC& fdc,
                       dealii::FullMatrix<double> &local_matrix)
        {
          assert(this->GetPart() == "New");
          this->GetProblem().BoundaryMatrix(fdc, local_matrix, 0., 1.);
        }
    private:
      InitialProblem<ForwardEulerProblem<OPTPROBLEM, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>, VECTOR, dealdim> * initial_problem_;
    };
}

#endif