crank_nicolson_problem.h

Go to the documentation of this file.

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

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

      std::string
      GetName()
      {
        return "Crank-Nicolson";
      }
      /******************************************************/

        InitialProblem<
            CrankNicolsonProblem<OPTPROBLEM, SPARSITYPATTERN, VECTOR, dopedim,
                dealdim, FE, DH>, VECTOR, dealdim>&
        GetInitialProblem()
        {
          if (initial_problem_ == NULL)
          {
            initial_problem_ = new InitialProblem<
                CrankNicolsonProblem<OPTPROBLEM, SPARSITYPATTERN, VECTOR,
                    dopedim, dealdim, FE, DH>, VECTOR, dealdim>(*this);
          }
          return *initial_problem_;
        }
    
      /******************************************************/
      CrankNicolsonProblem<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;
              // The remaining parts; e.g. for fluid problems: laplace, convection, etc.
              // Multiplication by 1/2 due to CN discretization
              this->GetProblem().ElementEquation(edc, tmp, 0.5 * 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")
            {
              dealii::Vector<double> tmp(local_vector);
              tmp = 0.0;

              // The explicit parts with old_time_values; e.g. for fluid problems: laplace, convection, etc.
              // Multiplication by 1/2 due to CN discretization
              this->GetProblem().ElementEquation(edc, tmp, 0.5 * 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")
            {
              this->GetProblem().ElementRhs(edc, local_vector, 0.5 * scale);
            }
          else if (this->GetPart() == "Old")
            {
              this->GetProblem().ElementRhs(edc, local_vector, 0.5 * 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")
           {
             this->GetProblem().PointRhs(param_values, domain_values, rhs_vector,
                 0.5 * scale);
           }
           else if (this->GetPart() == "Old")
           {
             this->GetProblem().PointRhs(param_values, domain_values, rhs_vector,
                 (0.5) * 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);

          // multiplication with 1/2 for scale due to CN discretization,
          //no multiplication with 1/2 for scale_ico due to implicit treatment of pressure, etc. (in the case of fluid problems)
          this->GetProblem().ElementMatrix(edc, local_matrix, 0.5, 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.5 * scale, scale);
            }
          else if (this->GetPart() == "Old")
            {
              this->GetProblem().FaceEquation(fdc, local_vector, 0.5 * 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.5 * scale, scale);
            }
          else if (this->GetPart() == "Old")
            {
              this->GetProblem().InterfaceEquation(fdc,  local_vector,
                                                   0.5 * 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");
          // Hier nicht mit this->GetProblem().GetBaseProblem().GetSpaceTimeHandler()->GetStepSize() multiplizieren, da local_matrix schon skaliert ist
          dealii::FullMatrix<double> m(local_matrix);

          m = 0.;
          // Multiplication with 1/2 due to CN time discretization
          this->GetProblem().FaceMatrix(fdc, m, 0.5, 1.);

          local_matrix.add(1., m);
        }

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

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

          m = 0.;
          // Multiplication with 1/2 due to CN time discretization
          this->GetProblem().InterfaceMatrix(fdc,  m, 0.5, 1.);

          local_matrix.add(1., m);
        }

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

     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.5 * scale, scale);
            }
          else if (this->GetPart() == "Old")
            {
              this->GetProblem().BoundaryEquation(fdc, local_vector, 0.5 * 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");
          dealii::FullMatrix<double> m(local_matrix);

          m = 0.;
          // Multiplication with 1/2 due to CN time discretization
          this->GetProblem().BoundaryMatrix(fdc, m, 0.5, 1.);
          local_matrix.add(1., m);

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

#endif