higher_order_dwrc_control.h

Go to the documentation of this file.

#ifndef HIGHER_ORDER_DWRC_CONTROL_H_
#define HIGHER_ORDER_DWRC_CONTROL_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 HigherOrderDWRContainerControl : public DWRDataContainer<STH, IDC, EDC, FDC,
        VECTOR>
    {
      public:
        HigherOrderDWRContainerControl(STH& higher_order_sth, IDC& higher_order_idc,
            DOpEtypes::VectorStorageType control_behavior, 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), PI_h_q_(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)
          {
            PI_h_q_ = new ControlVector<VECTOR>(&GetHigherOrderSTH(),
                                                control_behavior);
          }
        }

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

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

        template<class STH2>
          void
          Initialize(STH2* sth, unsigned int control_n_blocks,
                     std::vector<unsigned int>& control_block_component, 
                     unsigned int state_n_blocks,
                     std::vector<unsigned int>& state_block_component)
          {
            sth_ = dynamic_cast<STH*>(sth);
            control_n_blocks_        =  control_n_blocks;
            control_block_component_ = &control_block_component;
            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()
        {
          return *PI_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)
        {
          VECTOR q_high;
          q_high.reinit(GetPI_h_q().GetSpacialVector());

          dealii::FETools::extrapolate(
              GetSTH().GetControlDoFHandler().GetDEALDoFHandler(),
              q.GetSpacialVector(),
              GetHigherOrderSTH().GetControlDoFHandler().GetDEALDoFHandler(),
              GetHigherOrderSTH().GetControlDoFConstraints(),
              GetPI_h_q().GetSpacialVector());

          dealii::FETools::interpolate(
              GetSTH().GetControlDoFHandler().GetDEALDoFHandler(),
              q.GetSpacialVector(),
              GetHigherOrderSTH().GetControlDoFHandler().GetDEALDoFHandler(),
              GetHigherOrderSTH().GetControlDoFConstraints(), q_high);

          GetPI_h_q().GetSpacialVector().add(-1., q_high);
          //FIXME With this construction we can not deal with control contraints,
          //There, the real weight needs to be build on the elements...
        }

        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 control_n_blocks_, state_n_blocks_;
        std::vector<unsigned int>* control_block_component_;
        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_;
        ControlVector<VECTOR> * PI_h_q_;
    };

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

      GetHigherOrderSTH().ReInit(control_n_blocks_, *control_block_component_, state_n_blocks_, *state_block_component_);
      if (this->GetEETerms() == DOpEtypes::primal_only
          || this->GetEETerms() == DOpEtypes::mixed
          || this->GetEETerms() == DOpEtypes::mixed_control)
      {
        GetPI_h_z().ReInit();
      }
      if (this->GetEETerms() == DOpEtypes::dual_only
          || this->GetEETerms() == DOpEtypes::mixed
          || this->GetEETerms() == DOpEtypes::mixed_control)
      {
        GetPI_h_u().ReInit();
      }
      if( this->GetEETerms() == DOpEtypes::mixed_control)
      {
        GetPI_h_q().ReInit();
      }
    }
}

#endif /* HIGHER_ORDER_DWRC_H_ */