trueerrorestimator.hh

#ifndef __TRUEERROR_ESTIMATOR_HH__
#define __TRUEERROR_ESTIMATOR_HH__
 
// Local includes
#include "estimatorinterface.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class TrueSolutionFunction, class Norm>
    class TrueErrorEstimator
      : public DefaultEstimator<TrueErrorEstimator<TrueSolutionFunction, Norm> >
    {
      typedef TrueSolutionFunction TrueSolutionFunctionType;
      typedef Norm                 NormType;
      typedef TrueErrorEstimator<TrueSolutionFunctionType, NormType> ThisType;
      typedef DefaultEstimator<ThisType> BaseType;
     public:
      // Interface types
      typedef typename TrueSolutionFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::RangeFieldType RangeFieldType;
      typedef typename DiscreteFunctionSpaceType::GridPartType   GridPartType;
      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 DiscreteFunctionSpaceType::IteratorType GridIteratorType;
      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)
        : trueSolution_(trueSolution),
          dfSpace_(trueSolution_.space()),
          gridPart_(const_cast<GridPartType&>(dfSpace_.gridPart())),
          indexSet_(gridPart_.indexSet()),
          grid_(gridPart_.grid()),
          norm_(gridPart_)
      {
        localIndicators_.resize(indexSet_.size(0));
      }
 
      template<class DiscreteFunctionType>
      RangeFieldType estimate(const DiscreteFunctionType& uh)
      {
        // clear all local estimators
        clear();
 
        // calculate local estimator
        const GridIteratorType endGrid = dfSpace_.end();
        for (GridIteratorType it = dfSpace_.begin(); it != endGrid; ++it) {
          estimateLocal(*it, uh);
        }
 
        double error = 0.0;
 
        // sum up local estimators
        const IteratorType endEstimator = end();
        for (IteratorType it = begin(); it != endEstimator; ++it) {
          error += *it;
        }
 
        // obtain global sum
        error = grid_.comm().sum(error);
 
        // print error estimator
        std::cout << "Computed Error: " << std::sqrt(error) << std::endl;
        return std::sqrt(error);
      }
 
      const RangeFieldType& operator[](const size_t& idx) const {
        return localIndicators_[idx];
      }
 
      const RangeFieldType& operator[](const ElementType& entity) const {
        return (*this)[indexSet_.index(entity)];
      }
 
      size_t size() const {
        return indexSet_.size(0);
      }
 
      IteratorType begin() const {
        return localIndicators_.begin();
      }
 
      IteratorType end() const {
        return localIndicators_.end();
      }
 
      GridType& grid() const {
        return grid_;
      }
 
      const DiscreteFunctionSpaceType& space() const
      {
        return dfSpace_;
      }
 
     private:
      void clear ()
      {
        // resize and clear
        localIndicators_.resize(indexSet_.size(0));
        const typename ErrorIndicatorType::iterator end = localIndicators_.end();
        for (typename ErrorIndicatorType::iterator  it = localIndicators_.begin();
             it != end; ++it ) {
          *it = 0.0;
        }
      }
 
      template<class DiscreteFunctionType>
      void estimateLocal(const ElementType &entity, const DiscreteFunctionType& uh)
      {
        const int index = indexSet_.index(entity);
 
        // Crazy L2-norm takes a 1d vector as argument. Why?
        FieldVector<RangeFieldType, 1> dist(0);
        norm_.distanceLocal(entity, 2*uh.space().order()+2, trueSolution_.localFunction(entity), uh.localFunction(entity), dist);
        localIndicators_[index] = dist;
      }
     private:
      const TrueSolutionFunctionType& trueSolution_;
      const DiscreteFunctionSpaceType& dfSpace_;
      GridPartType& gridPart_;
      const IndexSetType& indexSet_;
      GridType& grid_;
      ErrorIndicatorType localIndicators_;
      NormType norm_;
    };
 
    template<class TrueSolutionFunction, class Norm>
    struct EstimatorTraits<TrueErrorEstimator<TrueSolutionFunction, Norm> >
    {
      typedef typename TrueSolutionFunction::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename
      std::vector<typename DiscreteFunctionSpaceType::RangeFieldType>::const_iterator
      IteratorType;
    };
 
 
 
  } // ACFem
 
} // Dune
 
#endif // __ELLIPTIC_ESTIMATOR_HH__