functionutilities.hh

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=8 sw=2 sts=2:
 
#ifndef DUNE_PDELAB_COMMON_FUNCTIONUTILITIES_HH
#define DUNE_PDELAB_COMMON_FUNCTIONUTILITIES_HH
 
#include <limits>
#include <ostream>
#include <memory>
 
#include <dune/common/shared_ptr.hh>
#include <dune/common/debugstream.hh>
 
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/type.hh>
 
#include <dune/grid/common/gridenums.hh>
#include <dune/grid/utility/hierarchicsearch.hh>
 
namespace Dune {
  namespace PDELab {
 
 
 
 
    template<typename GF>
    typename GF::Traits::RangeType integrateGridFunction(const GF& gf,
                               unsigned qorder = 1) {
      typename GF::Traits::RangeType sum;
      integrateGridFunction(gf, sum, qorder);
      return sum;
    }
 
 
    template<typename GF>
    void integrateGridFunction(const GF& gf,
                               typename GF::Traits::RangeType& sum,
                               unsigned qorder = 1) {
      typedef typename GF::Traits::GridViewType GV;
      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;
 
      sum = 0;
      Range val;
      for(const auto& cell : elements(gf.getGridView(),Dune::Partitions::interior)) {
        const Geometry& geo = cell.geometry();
        Dune::GeometryType gt = geo.type();
        const QR& rule = QRs::rule(gt,qorder);
        for (const auto& qip : rule) {
          // evaluate the given grid functions at integration point
          gf.evaluate(cell,qip.position(),val);
 
          // accumulate error
          val *= qip.weight() * geo.integrationElement(qip.position());
          sum += val;
        }
      }
    }
 
 
    template<typename GF>
    class GridFunctionProbe {
      typedef typename GF::Traits::GridViewType GV;
      typedef typename GV::template Codim<0>::Entity Entity;
      typedef typename GF::Traits::DomainType Domain;
      typedef typename GF::Traits::RangeType Range;
 
    public:
 
      GridFunctionProbe(const GF& gf, const Domain& xg)
        : gfp(Dune::stackobject_to_shared_ptr(gf))
      {
        hierarchic_search(gfp->getGridView(), xg);
      }
 
      GridFunctionProbe(std::shared_ptr<const GF> gf, const Domain& xg)
        : gfp(gf)
      {
        hierarchic_search(gfp->getGridView(), xg);
      }
 
      GridFunctionProbe(const GV& gv, const Domain& xg)
      {
        hierarchic_search(gv, xg);
      }
 
 
      void setGridFunction(const GF &gf) {
        gfp = Dune::stackobject_to_shared_ptr(gf);
      }
 
 
      void setGridFunction(const GF *gf) {
        gfp.reset(gf);
      }
 
 
      void setGridFunction(const std::shared_ptr<const GF> &gf) {
        gfp = gf;
      }
 
 
      void eval_all(Range& val) const {
        typedef typename GF::Traits::RangeFieldType RF;
        if(evalRank == gfp->getGridView().comm().size())
          val = std::numeric_limits<RF>::quiet_NaN();
        else {
          if(gfp->getGridView().comm().rank() == evalRank)
            gfp->evaluate(*e, xl, val);
          gfp->getGridView().comm().broadcast(&val,1,evalRank);
        }
      }
 
 
      void eval(Range& val, int rank = 0) const {
        eval_all(val);
      }
 
 
      int getEvalRank() {
        return evalRank;
      }
 
    private:
      void hierarchic_search(const GV& gv, const Domain& xg)
      {
        xl = 0;
        evalRank = gv.comm().size();
        int myRank = gv.comm().rank();
        try {
          e.reset(new Entity
                  (HierarchicSearch<typename GV::Grid, typename GV::IndexSet>
                   (gv.grid(), gv.indexSet()).
                   template findEntity<Interior_Partition>(xg)));
          // make sure only interior entities are accepted
          if(e->partitionType() == InteriorEntity)
            evalRank = myRank;
        }
        catch(const Dune::GridError&) { /* do nothing */ }
        evalRank = gv.comm().min(evalRank);
        if(myRank == evalRank)
          xl = e->geometry().local(xg);
        else
          e.reset();
        if(myRank == 0 && evalRank == gv.comm().size())
          dwarn << "Warning: GridFunctionProbe at (" << xg << ") is outside "
                << "the grid" << std::endl;
      }
 
      std::shared_ptr<const GF> gfp;
      std::shared_ptr<Entity> e;
      Domain xl;
      int evalRank;
    };
 
 
  } // namespace PDELab
} // namespace Dune
 
#endif // DUNE_PDELAB_COMMON_FUNCTIONUTILITIES_HH