functionutilities.hh

// -*- tab-width: 2; indent-tabs-mode: nil -*-
#ifndef DUNE_NPDE_FUNCTIONUTILITIES_HH
#define DUNE_NPDE_FUNCTIONUTILITIES_HH

#include<dune/geometry/quadraturerules.hh>
#include<dune/geometry/type.hh>

#include <dune/grid/common/gridenums.hh>
#include <dune/grid/utility/hierarchicsearch.hh>

template<typename GF>
void getGridFunctionMaxMin(const GF& gf,
                           typename GF::Traits::RangeType& maxx,
                           typename GF::Traits::RangeType& minn,
                           unsigned qorder = 1) {
  typedef typename GF::Traits::GridViewType GV;
  typedef typename GV::template Codim<0>::
    template Partition<Dune::Interior_Partition>::Iterator EIterator;
  typedef typename GV::template Codim<0>::Geometry Geometry;
  typedef typename GF::Traits::RangeType Range;
  typedef typename GF::Traits::DomainFieldType DF;
  static const int dimD = GF::Traits::dimDomain;
  typedef Dune::QuadratureRule<DF,dimD> QR;
  typedef Dune::QuadratureRules<DF,dimD> QRs;
  typedef typename QR::const_iterator QIterator;

  maxx = 0; minn= 0;
  Range val;
  const EIterator eend = gf.getGridView().template end
    <0,Dune::Interior_Partition>();
  for(EIterator eit = gf.getGridView().template begin<0,
        Dune::Interior_Partition>(); eit != eend; ++eit) {
    const Geometry& geo = eit->geometry();
    Dune::GeometryType gt = geo.type();
    const QR& rule = QRs::rule(gt,qorder);
    const QIterator qend = rule.end();

    for (QIterator qit=rule.begin(); qit != qend; ++qit)
      {
        // evaluate the given grid functions at integration point
        gf.evaluate(*eit,qit->position(),val);

        maxx = std::max(val[0],maxx[0]);
        minn = std::min(val[0],minn[0]);
      }
  }
}


/*! \brief Adapter returning f1(x)-f2(x) for two given grid functions

  \tparam T1  a grid function type
  \tparam T2  a grid function type
*/
template<typename T1, typename T2>
class DifferenceAdapter
  : public Dune::PDELab::GridFunctionBase<
  Dune::PDELab::GridFunctionTraits<typename T1::Traits::GridViewType,
                                   typename T1::Traits::RangeFieldType,
                                   1,Dune::FieldVector<typename T1::Traits::RangeFieldType,1> >
  ,DifferenceAdapter<T1,T2> >
{
public:
  typedef Dune::PDELab::GridFunctionTraits<typename T1::Traits::GridViewType,
                                           typename T1::Traits::RangeFieldType,
                                           1,Dune::FieldVector<typename T1::Traits::RangeFieldType,1> > Traits;

  //! constructor
  DifferenceAdapter (const T1& t1_, const T2& t2_) : t1(t1_), t2(t2_) {}