errorest.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2013-2015 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef ERROREST_HH
#define ERROREST_HH
 
#include <vector>
 
#include <boost/fusion/include/for_each.hpp>
#include <boost/utility.hpp>
#include <boost/multi_array.hpp>
 
#include "grid/common/rangegenerators.hh"
 
namespace Kaskade 
{
  class RefinementCriterion
  {
  public:
    std::vector<double> threshold(boost::multi_array<double,2> const& normalizedErrors) const;
    
  private:
    virtual double computeThreshold(std::vector<double>& normalizedErrors, double totalError, int j) const = 0;
  };
  
  template <class GridManager, class GridView>
  void markAndRefine(GridManager& gridManager, GridView const& gridView, boost::multi_array<double,2> const& normalizedErrors, 
                     std::vector<double> const& threshold)
  {
    for (auto const& cell: Dune::elements(gridView))
    {
      auto idx = gridView.indexSet().index(cell);
      for (int j=0; j<threshold.size(); ++j)
        if (normalizedErrors[idx][j] >= threshold[j])
        {
          gridManager.mark(1,cell);
          break;
        }
    }
    gridManager.adaptAtOnce();
  }
  
  class FixedFractionCriterion: public RefinementCriterion
  {
  public:
    FixedFractionCriterion(double fraction = 0.2);
    
  private:
    double fraction;
    virtual double computeThreshold(std::vector<double>& normalizedErrors, double totalError, int j) const;
  };
  
  class BulkCriterion: public RefinementCriterion
  {
  public:
    BulkCriterion(double fraction = 0.2);
    
  private:
    double fraction;
    virtual double computeThreshold(std::vector<double>& normalizedErrors, double totalError, int j) const;
  };
  
  class MaxValueCriterion: public RefinementCriterion
  {
  public:
    MaxValueCriterion(double fraction = 0.2);
    
  private:
    double fraction;
    virtual double computeThreshold(std::vector<double>& normalizedErrors, double totalError, int j) const;
  };
  
  class BabuskaRheinboldtCriterion: public MaxValueCriterion
  {
  public:
    BabuskaRheinboldtCriterion(std::vector<int> const& order);
    
  private:
    std::vector<int> order;
    virtual double computeThreshold(std::vector<double>& normalizedErrors, double totalError, int j) const;
  };
  
  //---------------------------------------------------------------------------------------------
  
  namespace ErrorestDetail
  {
    template <class GroupByCell>
    struct GroupedSummationCollector
    {
      GroupedSummationCollector(GroupByCell const& group_):
        group(group_)
      {}
      
      template <class Cell>
      int integrationOrder(Cell const& , int shapeFunctionOrder) const
      {
        return shapeFunctionOrder;
      }
      
      // need to define this because boost::multi_array does not perform assignment if the arrays have
      // different shape.
      GroupedSummationCollector<GroupByCell>& operator=(GroupedSummationCollector<GroupByCell> const& c)
      {
        auto shape = c.sums.shape();
 
        sums.resize(boost::extents[shape[0]][shape[1]]);
        sums = c.sums;
        group = c.group;
 
        return *this;
      }
      
      template <class Cell, class Index, class LocalPosition, class Sequence>
      void operator()(Cell const& /* cell */, Index idx,
                      LocalPosition const& /* xi */, double weight,
                      Sequence const& x)
      {
        if (sums.num_elements()==0)
          sums.resize(boost::extents[group.nGroups][boost::fusion::size(x)]);
        int i = 0; // for_each needs constant functor...
        boost::fusion::for_each(x,Add(sums,group[idx],i,weight));
      }
      
      void join(GroupedSummationCollector<GroupByCell> const& c)
      {
        auto shape = c.sums.shape();
        auto myshape = sums.shape();
        
        if(c.sums.num_elements()==0)  // no data in collector c
          return;
 
        if (sums.num_elements()==0 || myshape[0]!=shape[0] || myshape[1]!=shape[1]) // I'm not yet initialized
        {
          sums.resize(boost::extents[shape[0]][shape[1]]);
          sums = c.sums;                                  // -> do a simple copy
        }
        else
          for (size_t i=0; i<shape[0]; ++i)               // otherwise add up
            for (size_t j=0; j<shape[1]; ++j)
              sums[i][j] += c.sums[i][j];
      }
 
      boost::multi_array<double,2> sums;
 
    private:
      GroupByCell group;
 
      struct Add
      {
        Add(boost::multi_array<double,2>& sums_, int idx_, int& i_, double w_):
          sums(sums_), idx(idx_), i(i_) , w(w_)
        {}
 
        template <class T>
        void operator()(T const& t) const { sums[idx][i++] += w*t; }
 
      private:
        boost::multi_array<double,2>& sums;
        int idx;
        int& i;
        double w;
      };
    };
    
  }
}
 
#endif