trueerrorestimator.hh

#ifndef __TRUEERROR_ESTIMATOR_HH__
#define __TRUEERROR_ESTIMATOR_HH__
 
#include <numeric>
 
#include <dune/fem/function/localfunction/const.hh>
 
#include "../common/entitystorage.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class TrueSolutionFunction, class Norm>
    class TrueErrorEstimator
      : public EntityStorage<typename TrueSolutionFunction::GridPartType,
                             typename TrueSolutionFunction::GridPartType::ctype>
    {
      typedef TrueSolutionFunction TrueSolutionFunctionType;
      typedef Norm                 NormType;
      using ThisType = TrueErrorEstimator<TrueSolutionFunctionType, NormType>;
      using BaseType = EntityStorage<typename TrueSolutionFunction::GridPartType,
                                     typename TrueSolutionFunction::GridPartType::ctype>;
     public:
      // Interface types
      typedef typename TrueSolutionFunctionType::GridPartType GridPartType;
      typedef typename TrueSolutionFunction::FunctionSpaceType::RangeFieldType RangeFieldType;
      typedef typename GridPartType::GridType              GridType;
      typedef typename GridType::template Codim<0>::Entity ElementType;
     protected:
      // Could live without that, but we cache the results of the local
      // error computations in a vector.
      typedef std::vector<RangeFieldType> ErrorIndicatorType;
      typedef typename GridPartType::IndexSetType IndexSetType;
     public:
      typedef typename ErrorIndicatorType::const_iterator IteratorType;
 
      // Constructor. The estimator-interface expects the estimate-method
      // to accept only one argument -- the discrete solution to be
      // estimated -- so the "true" function has to be passed on through
      // the estimator (or: at least this is the easiest way)
      explicit TrueErrorEstimator(const TrueSolutionFunctionType& trueSolution)
        : BaseType(trueSolution.gridPart()),
          trueSolution_(trueSolution),
          gridPart_(trueSolution.gridPart()),
          indexSet_(gridPart_.indexSet()),
          grid_(gridPart_.grid()),
          norm_(gridPart_)
      {
      }
 
      template<class DiscreteFunctionType>
      RangeFieldType estimate(const DiscreteFunctionType& uh)
      {
        // clear all local estimators
        clear();
 
        Fem::ConstLocalFunction<DiscreteFunctionType> uhLocal(uh);
 
        // calculate local estimator
        const auto endGrid = uh.space().end();
        for (auto it = uh.space().begin(); it != endGrid; ++it) {
          estimateLocal(*it, uhLocal);
        }
 
        RangeFieldType error = 0.0;
 
        // sum up local estimators
        error = std::accumulate(this->storage_.begin(), this->storage_.end(), error);
 
        // obtain global sum
        error = grid_.comm().sum(error);
 
        // print error estimator
        std::cout << "Computed Error: " << std::sqrt(error) << std::endl;
        return std::sqrt(error);
      }
 
     private:
      void clear ()
      {
        // resize and clear
        BaseType::clear();
      }
 
      template<class DiscreteLocalFunctionType>
      void estimateLocal(const ElementType &entity, DiscreteLocalFunctionType& uh)
      {
        const int index = indexSet_.index(entity);
 
        uh.bind(entity);
        trueSolution_.bind(entity);
 
        // Crazy L2-norm takes a 1d vector as argument. Why?
        FieldVector<RangeFieldType, 1> dist(0);
        norm_.distanceLocal(entity, 2*uh.order()+2, trueSolution_, uh, dist);
        this->storage_[index] = dist;
 
        trueSolution_.unbind();
        uh.unbind();
      }
     private:
      Fem::ConstLocalFunction<TrueSolutionFunctionType> trueSolution_;
      const GridPartType& gridPart_;
      const IndexSetType& indexSet_;
      const GridType& grid_;
      NormType norm_;
    };
 
 
 
  } // ACFem
 
} // Dune
 
#endif // __ELLIPTIC_ESTIMATOR_HH__