backward_euler_problem.h

Go to the documentation of this file.

#ifndef BackwardEulerProblem_H_
#define BackwardEulerProblem_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 BackwardEulerProblem : public PrimalTSBase<OPTPROBLEM,
        SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>
    {
      public:
        BackwardEulerProblem(OPTPROBLEM& OP) :
            PrimalTSBase<OPTPROBLEM, SPARSITYPATTERN, VECTOR, dopedim, dealdim,
                FE, DH>(OP)
        {
          initial_problem_ = NULL;
        }
        ~BackwardEulerProblem()
        {
          if (initial_problem_ != NULL)
            delete initial_problem_;
        }

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

        std::string
        GetName()
        {
          return "backward Euler";
        }
        /******************************************************/
      
        InitialProblem<
            BackwardEulerProblem<OPTPROBLEM, SPARSITYPATTERN, VECTOR, dopedim,
                dealdim, FE, DH>, VECTOR, dealdim>&
        GetInitialProblem()
        {
          if (initial_problem_ == NULL)
          {
            initial_problem_ = new InitialProblem<
                BackwardEulerProblem<OPTPROBLEM, SPARSITYPATTERN, VECTOR,
                    dopedim, dealdim, FE, DH>, VECTOR, dealdim>(*this);
          }
          return *initial_problem_;
        }

        /******************************************************/
        
        BackwardEulerProblem<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,
                                                 scale,
                                                 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")
            {
              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")
            {
              this->GetProblem().ElementRhs(edc, local_vector,
                                            scale);
            }
            else if (this->GetPart() == "Old")
            {
            }
            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")
          {
            this->GetProblem().PointRhs(param_values, domain_values, rhs_vector,
                                        scale);
          }
          else if (this->GetPart() == "Old")
          {
          }
          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,
                                             1.,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,
                                              scale,
                                              scale);
            }
            else if (this->GetPart() == "Old")
            {
            }
            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,
                                                   scale,
                                                   scale);
            }
            else if (this->GetPart() == "Old")
            {
            }
            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,
                                          1.,1.);

          }

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

         template<typename FDC>
          void
          InterfaceMatrix(const FDC& fdc,
              dealii::FullMatrix<double> &local_matrix)
          {
            assert(this->GetPart() == "New");
            this->GetProblem().InterfaceMatrix(fdc, local_matrix,
                                               1.,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,
                                                  scale, scale);
            }
            else if (this->GetPart() == "Old")
            {
            }
            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,
                                              1., 1.);
          }

      private:
        InitialProblem<
            BackwardEulerProblem<OPTPROBLEM, SPARSITYPATTERN, VECTOR, dopedim,
                dealdim, FE, DH>, VECTOR, dealdim> * initial_problem_;
    };
}

#endif