higher_order_dwrc.h

Go to the documentation of this file.

#ifndef HIGHER_ORDER_DWRC_H_
#define HIGHER_ORDER_DWRC_H_

#include "dwrdatacontainer.h"
#include <deal.II/fe/fe_tools.h>

namespace DOpE
{
  template<class STH, class IDC, class EDC, class FDC, typename VECTOR>
    class HigherOrderDWRContainer : public DWRDataContainer<STH, IDC, EDC, FDC,
        VECTOR>
    {
      public:
        HigherOrderDWRContainer(STH& higher_order_sth, IDC& higher_order_idc,
            DOpEtypes::VectorStorageType state_behavior, ParameterReader &param_reader,
            DOpEtypes::EETerms ee_terms = DOpEtypes::EETerms::mixed,
            DOpEtypes::ResidualEvaluation res_eval = DOpEtypes::strong_residual)
            : DWRDataContainer<STH, IDC, EDC, FDC, VECTOR>(ee_terms), sth_higher_order_(
              higher_order_sth), idc_higher_order_(higher_order_idc), res_eval_(res_eval),
              PI_h_u_(NULL), PI_h_z_(NULL)
        {
          if (this->GetEETerms() == DOpEtypes::primal_only
              || this->GetEETerms() == DOpEtypes::mixed
              || this->GetEETerms() == DOpEtypes::mixed_control)
          {
            PI_h_z_ = new StateVector<VECTOR>(&GetHigherOrderSTH(),
                state_behavior, param_reader);
          }
          if (this->GetEETerms() == DOpEtypes::dual_only
              || this->GetEETerms() == DOpEtypes::mixed
              || this->GetEETerms() == DOpEtypes::mixed_control)
          {
            PI_h_u_ = new StateVector<VECTOR>(&GetHigherOrderSTH(),
                state_behavior, param_reader);
          }
          if (this->GetEETerms() == DOpEtypes::mixed_control)
          {
            throw DOpEException("Wrong Constructor for type 'mixed_control'",
                  "HigherOrderDWRContainer::HigherOrderDWRContainer");
          }
        }

        virtual
        ~HigherOrderDWRContainer()
        {
          if (PI_h_z_ != NULL)
            delete PI_h_z_;
          if (PI_h_u_ != NULL)
            delete PI_h_u_;
        }

        std::string
        GetName() const
        {
          return "DWR-Estimator";
        }

        template<class STH2>
          void
          Initialize(STH2* sth, unsigned int state_n_blocks,
              std::vector<unsigned int>& state_block_component)
          {
            _sth = dynamic_cast<STH*>(sth);
            _state_n_blocks = state_n_blocks;
            _state_block_component = &state_block_component;
          }

        void
        ReInit(unsigned int n_elements);

        STH&
        GetWeightSTH()
        {
          return GetHigherOrderSTH();
        }

        const STH&
        GetWeightSTH() const
        {
          return GetHigherOrderSTH();
        }

        IDC&
        GetWeightIDC()
        {
          return GetHigherOrderIDC();
        }

        const IDC&
        GetWeightIDC() const
        {
          return GetHigherOrderIDC();
        }

        StateVector<VECTOR>&
        GetPI_h_u()
        {
          return *PI_h_u_;
        }

        StateVector<VECTOR>&
        GetPI_h_z()
        {
          return *PI_h_z_;
        }
        ControlVector<VECTOR>&
        GetPI_h_q()
        {
          throw DOpEException("There is no Control in PDE Problems!",
            "HigherOrderDWRContainer::GetPI_h_q");
        }

        void
        PreparePI_h_u(const StateVector<VECTOR>& u)
        {
          VECTOR u_high;
          u_high.reinit(GetPI_h_u().GetSpacialVector());

          dealii::FETools::extrapolate(
              GetSTH().GetStateDoFHandler().GetDEALDoFHandler(),
              u.GetSpacialVector(),
              GetHigherOrderSTH().GetStateDoFHandler().GetDEALDoFHandler(),
              GetHigherOrderSTH().GetStateDoFConstraints(),
              GetPI_h_u().GetSpacialVector());
          dealii::FETools::interpolate(
              GetSTH().GetStateDoFHandler().GetDEALDoFHandler(),
              u.GetSpacialVector(),
              GetHigherOrderSTH().GetStateDoFHandler().GetDEALDoFHandler(),
              GetHigherOrderSTH().GetStateDoFConstraints(), u_high);
          GetPI_h_u().GetSpacialVector().add(-1., u_high);
        }

        void
        PreparePI_h_u(const VECTOR& u)
        {
          VECTOR u_high;
          u_high.reinit(GetPI_h_u().GetSpacialVector());

          dealii::FETools::extrapolate(
              GetSTH().GetStateDoFHandler().GetDEALDoFHandler(), u,
              GetHigherOrderSTH().GetStateDoFHandler().GetDEALDoFHandler(),
              GetHigherOrderSTH().GetStateDoFConstraints(),
              GetPI_h_u().GetSpacialVector());
          dealii::FETools::interpolate(
              GetSTH().GetStateDoFHandler().GetDEALDoFHandler(), u,
              GetHigherOrderSTH().GetStateDoFHandler().GetDEALDoFHandler(),
              GetHigherOrderSTH().GetStateDoFConstraints(), u_high);
          GetPI_h_u().GetSpacialVector().add(-1., u_high);
        }

        void
        PreparePI_h_z(const StateVector<VECTOR>& z)
        {
          VECTOR z_high;
          z_high.reinit(GetPI_h_z().GetSpacialVector());

          dealii::FETools::extrapolate(
              GetSTH().GetStateDoFHandler().GetDEALDoFHandler(),
              z.GetSpacialVector(),
              GetHigherOrderSTH().GetStateDoFHandler().GetDEALDoFHandler(),
              GetHigherOrderSTH().GetStateDoFConstraints(),
              GetPI_h_z().GetSpacialVector());

          dealii::FETools::interpolate(
              GetSTH().GetStateDoFHandler().GetDEALDoFHandler(),
              z.GetSpacialVector(),
              GetHigherOrderSTH().GetStateDoFHandler().GetDEALDoFHandler(),
              GetHigherOrderSTH().GetStateDoFConstraints(), z_high);

          GetPI_h_z().GetSpacialVector().add(-1., z_high);

        }
        void
          PreparePI_h_q(const ControlVector<VECTOR>& /*q*/)
        {
          throw DOpEException("There is no Control in PDE Problems!",
            "HigherOrderDWRContainer::PreparePI_h_q");

        }

        virtual EDC&
        GetElementWeight() const
        {
          return GetHigherOrderIDC().GetElementDataContainer();
        }

        virtual FDC&
        GetFaceWeight() const
        {
          return GetHigherOrderIDC().GetFaceDataContainer();
        }

        bool
        NeedDual() const
        {
          return true;
        }

        virtual DOpEtypes::WeightComputation
        GetWeightComputation() const
        {
          return DOpEtypes::higher_order_interpolation;
        }

        virtual DOpEtypes::ResidualEvaluation
        GetResidualEvaluation() const
        {
//          return DOpEtypes::strong_residual;
          return res_eval_;
        }

        inline void
          ResidualModifier(double& /*res*/)
        {
          
        }
        inline void
          VectorResidualModifier(dealii::Vector<double>& /*res*/)
        {
      
        }

      protected:
        STH&
        GetSTH()
        {
          return *_sth;
        }

        STH&
        GetHigherOrderSTH()
        {
          return sth_higher_order_;
        }

        const STH&
        GetHigherOrderSTH() const
        {
          return sth_higher_order_;
        }

        IDC&
        GetHigherOrderIDC()
        {
          return idc_higher_order_;
        }

        const IDC&
        GetHigherOrderIDC() const
        {
          return idc_higher_order_;
        }

      private:
        unsigned int _state_n_blocks;
        std::vector<unsigned int>* _state_block_component;

        STH& sth_higher_order_;
        STH* _sth;
        IDC& idc_higher_order_;

        const DOpEtypes::ResidualEvaluation res_eval_;

        StateVector<VECTOR> * PI_h_u_, *PI_h_z_;
    };

  template<class STH, class IDC, class EDC, class FDC, typename VECTOR>
    void
    HigherOrderDWRContainer<STH, IDC, EDC, FDC, VECTOR>::ReInit(
        unsigned int n_elements)
    {
      DWRDataContainer<STH, IDC, EDC, FDC, VECTOR>::ReInit(n_elements);

      GetHigherOrderSTH().ReInit(_state_n_blocks, *_state_block_component);
      if (this->GetEETerms() == DOpEtypes::primal_only
          || this->GetEETerms() == DOpEtypes::mixed)
      {
        GetPI_h_z().ReInit();
      }
      if (this->GetEETerms() == DOpEtypes::dual_only
          || this->GetEETerms() == DOpEtypes::mixed)
      {
        GetPI_h_u().ReInit();
      }
    }
}

#endif /* HIGHER_ORDER_DWRC_H_ */