residualestimator.h

Go to the documentation of this file.

#ifndef RESIDUAL_ERROR_H_
#define RESIDUAL_ERROR_H_

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

namespace DOpE
{
  template<typename VECTOR>
    class ResidualErrorContainer : public DWRDataContainerBase<VECTOR>
    {
      public:
        ResidualErrorContainer(DOpEtypes::EETerms ee_terms =
            DOpEtypes::EETerms::mixed) :
            DWRDataContainerBase<VECTOR>(ee_terms)
        {
        }
        virtual void
        InitFace(double /*h*/) = 0;
        virtual void
        InitElement(double /*h*/) = 0;

    };

  template<class EDC, typename VECTOR>
    EDC*
    ExtractEDC(const ResidualErrorContainer<VECTOR>& /*dwrc*/)
    {
      return NULL;
    }
  template<class FDC, typename VECTOR>
    FDC*
    ExtractFDC(const ResidualErrorContainer<VECTOR>& /*dwrc*/)
    {
      return NULL;
    }
  template<class STH, typename VECTOR, int dim>
    class L2ResidualErrorContainer : public ResidualErrorContainer<VECTOR>
    {
      public:
        L2ResidualErrorContainer(STH& sth, DOpEtypes::VectorStorageType state_behavior,
            ParameterReader &param_reader, DOpEtypes::EETerms ee_terms =
                DOpEtypes::EETerms::mixed) :
            ResidualErrorContainer<VECTOR>(ee_terms), sth_(sth), PI_h_u_(NULL), PI_h_z_(
                NULL)
        {
          if (this->GetEETerms() == DOpEtypes::primal_only
              || this->GetEETerms() == DOpEtypes::mixed)
          {
            PI_h_z_ = new StateVector<VECTOR>(&GetSTH(), state_behavior,
                param_reader);
          }
          if (this->GetEETerms() == DOpEtypes::dual_only
              || this->GetEETerms() == DOpEtypes::mixed)
          {
            PI_h_u_ = new StateVector<VECTOR>(&GetSTH(), state_behavior,
                param_reader);
          }
          weight_ = 0.;
        }

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

        }

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

        void
        Initialize(unsigned int state_n_blocks,
            std::vector<unsigned int>& state_block_component)
        {
          state_n_blocks_ = state_n_blocks;
          state_block_component_ = &state_block_component;
        }

        void
        ReInit(unsigned int n_elements);

        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!",
              "L2ResidualErrorContainer::PreparePI_h_q");
        }

        void
        PreparePI_h_u(const StateVector<VECTOR>& /*u*/)
        {
          BuildConstantWeight(&(GetSTH().GetStateDoFHandler()),
              GetPI_h_u().GetSpacialVector());
        }

        void
        PreparePI_h_z(const StateVector<VECTOR>& /*z*/)
        {
          BuildConstantWeight(&(GetSTH().GetStateDoFHandler()),
              GetPI_h_z().GetSpacialVector());
        }
        void
        PreparePI_h_q(const ControlVector<VECTOR>& /*q*/)
        {
          throw DOpEException("There is no Control in PDE Problems!",
              "HigherOrderDWRContainer::PreparePI_h_q");

        }

        bool
        NeedDual() const
        {
          return false;
        }

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

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

        inline void
        ResidualModifier(double& res)
        {
          res = res * res * weight_;
        }

        inline void
        VectorResidualModifier(dealii::Vector<double>& res)
        {
          for (unsigned int i = 0; i < res.size(); i++)
            res(i) = res(i) * res(i) * weight_;
        }

        void
        InitFace(double h)
        {
          weight_ = h * h * h;
        }
        void
        InitElement(double h)
        {
          weight_ = h * h * h * h;
        }

      protected:
        STH&
        GetSTH()
        {
          return sth_;
        }

        template<template<int, int> class DH>
          void
          BuildConstantWeight(const DOpEWrapper::DoFHandler<dim, DH>* dofh,
              VECTOR& vals)
          {
            VectorTools::interpolate(sth_.GetMapping(),
                *(static_cast<const DH<dim, dim>*>(dofh)),
                ConstantFunction<dim>(1., dofh->get_fe().n_components()), vals);
          }

      private:
        unsigned int state_n_blocks_;
        std::vector<unsigned int>* state_block_component_;
        double weight_;

        STH& sth_;

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

  template<class STH, typename VECTOR, int dim>
    void
    L2ResidualErrorContainer<STH, VECTOR, dim>::ReInit(unsigned int n_elements)
    {
      DWRDataContainerBase<VECTOR>::ReInit(n_elements);

      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();
      }
    }

  /**************************************************************************/
  template<class STH, typename VECTOR, int dim>
    class H1ResidualErrorContainer : public ResidualErrorContainer<VECTOR>
    {
      public:
        H1ResidualErrorContainer(STH& sth, DOpEtypes::VectorStorageType state_behavior,
            ParameterReader &param_reader, DOpEtypes::EETerms ee_terms =
                DOpEtypes::EETerms::mixed) :
            ResidualErrorContainer<VECTOR>(ee_terms), sth_(sth), PI_h_u_(NULL), PI_h_z_(
                NULL)
        {
          if (this->GetEETerms() == DOpEtypes::primal_only
              || this->GetEETerms() == DOpEtypes::mixed)
          {
            PI_h_z_ = new StateVector<VECTOR>(&GetSTH(), state_behavior,
                param_reader);
          }
          if (this->GetEETerms() == DOpEtypes::dual_only
              || this->GetEETerms() == DOpEtypes::mixed)
          {
            PI_h_u_ = new StateVector<VECTOR>(&GetSTH(), state_behavior,
                param_reader);
          }
        }

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

        }

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

        void
        Initialize(unsigned int state_n_blocks,
            std::vector<unsigned int>& state_block_component)
        {
          state_n_blocks_ = state_n_blocks;
          state_block_component_ = &state_block_component;
        }

        void
        ReInit(unsigned int n_elements);

        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!",
              "H1ResidualErrorContainer::PreparePI_h_q");
        }

        void
        PreparePI_h_u(const StateVector<VECTOR>& /*u*/)
        {
          BuildConstantWeight(&(GetSTH().GetStateDoFHandler()),
              GetPI_h_u().GetSpacialVector());
        }

        void
        PreparePI_h_z(const StateVector<VECTOR>& /*z*/)
        {
          BuildConstantWeight(&(GetSTH().GetStateDoFHandler()),
              GetPI_h_z().GetSpacialVector());
        }
        void
        PreparePI_h_q(const ControlVector<VECTOR>& /*q*/)
        {
          throw DOpEException("There is no Control in PDE Problems!",
              "HigherOrderDWRContainer::PreparePI_h_q");

        }

        bool
        NeedDual() const
        {
          return false;
        }

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

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

        inline void
        ResidualModifier(double& res)
        {
          res = res * res * weight_;
        }
        inline void
        VectorResidualModifier(dealii::Vector<double>& res)
        {
          for (unsigned int i = 0; i < res.size(); i++)
            res(i) = res(i) * res(i) * weight_;
        }

        void
        InitFace(double h)
        {
          weight_ = h;
        }
        void
        InitElement(double h)
        {
          weight_ = h * h;
        }

      protected:
        STH&
        GetSTH()
        {
          return sth_;
        }

        template<template<int, int> class DH>
          void
          BuildConstantWeight(const DOpEWrapper::DoFHandler<dim, DH>* dofh,
              VECTOR& vals)
          {
            VectorTools::interpolate(sth_.GetMapping(),
                dofh->GetDEALDoFHandler(),
                ConstantFunction<dim>(1., dofh->get_fe().n_components()), vals);
          }

      private:
        unsigned int state_n_blocks_;
        std::vector<unsigned int>* state_block_component_;
        double weight_;

        STH& sth_;

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

  template<class STH, typename VECTOR, int dim>
    void
    H1ResidualErrorContainer<STH, VECTOR, dim>::ReInit(unsigned int n_elements)
    {
      DWRDataContainerBase<VECTOR>::ReInit(n_elements);

      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 /* RESIDUAL_ERROR */