integration.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) 2002-2019 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef FEM_INTEGRATION_HH
#define FEM_INTEGRATION_HH
 
#include <dune/geometry/quadraturerules.hh>
 
#include "fem/functionspace.hh"
#include "fem/forEach.hh"
 
namespace Kaskade
{
  namespace Integration_Detail
  {
    // C++ 14 has polymorphic lambdas..., but without method "order".
    template <class Function, class Functor>
    struct G
    {
      G(Function const& f_, Functor& g_)
      : f(f_), g(g_) {}
 
      template <class Cell, class QuadPoint, class Evaluator>
      void operator()(Cell const& cell, QuadPoint const& qp, Evaluator const& eval) const
      {
        g(cell,qp,f.value(eval));
      }
 
      template <class Evaluator>
      int order(Evaluator const& eval) const
      {
        return f.order(eval);
      }
 
      Function const& f;
      Functor& g;
    };
 
  }
 
  template <class Function, class Functor>
  void forEachQuadPoint(Function const& f, Functor& accumulator)
  {
    Integration_Detail::G<Function,Functor> gWrapper(f,accumulator);
  // this lambda would be great, but it doesn't provide method order
  //     auto gWrapper = [&] (auto const& cellPointer, auto const& quadPoint, auto w, auto const& eval)
  //     {
  //       accumulator(cellPointer,quadPoint,f.value(eval));
  //     };
    forEachQP(gWrapper,f.space());
  }
 
  template <class Functor, class Space>
  void forEachQP(Functor& g, Space const& space)
  {
    typedef typename Space::Grid Grid;
    typedef typename Grid::LeafGridView GridView;
    typedef typename Space::Evaluator Evaluator;
 
    ShapeFunctionCache<Grid,typename Space::Scalar> sfCache;
 
    GridView const& gridView = space.grid().leafGridView();
 
    Evaluator evaluator(space,&sfCache);
    for(auto const& cell : elements(gridView))
    {
      evaluator.moveTo(cell);
      Dune::GeometryType gt = cell.type();
      auto const& qr = Dune::QuadratureRules<typename Grid::ctype,Grid::dimension>::rule(gt,g.order(evaluator));
      evaluator.useQuadratureRule(qr,0);
 
      size_t nQuadPos = qr.size();
      for (size_t i=0; i<nQuadPos; ++i)
      {
        evaluator.evaluateAt(qr[i].position(),qr,i,0);
        g(cell,qr[i],evaluator);
      }
    }
  }
 
  template <class Functor, class Function, class Space>
  void forEachQuadPoint(Function const& f, Functor& g, Space& space)
  {
    typedef typename Space::Grid Grid;
    typedef typename Grid::LeafGridView GridView;
    typedef typename Space::Evaluator Evaluator;
 
    typedef ShapeFunctionCache<Grid,typename Space::Scalar> SfCache;
    SfCache sfCache;
 
    //   typedef typename IndexSet::template Codim<0>::template Partition<Dune::All_Partition> Entity;
    //   typedef typename Entity::Iterator CellIterator;
    //   IndexSet const& indexSet = space.grid().leafIndexSet();
    typedef typename GridView::template Codim<0>::Iterator CellIterator ;
    GridView const& gridView = f.space().grid().leafGridView() ;
 
    Evaluator evaluator(space, &sfCache);
 
    for (CellIterator ci=gridView.template begin<0,Dune::All_Partition>();
        ci!=gridView.template end<0,Dune::All_Partition>(); ++ci)
    {
      evaluator.moveTo(*ci);
 
      Dune::GeometryType gt = ci->type();
      auto const& qr = Dune::QuadratureRules<typename Grid::ctype,Grid::dimension>::rule(gt,f.order(evaluator));
 
      size_t nQuadPos = qr.size();
      for (size_t i=0; i<nQuadPos; ++i)
        g(typename GridView::template Codim<0>::EntityPointer(ci),qr[i],f.value(*ci,qr[i].position()));
    }
  }
 
  template <class FEFunction>
  typename FEFunction::ValueType integrateOverIntersection(FEFunction const& function,
                                                           typename FEFunction::Space::GridView::Intersection const& intersection,
                                                           typename FEFunction::Space::Evaluator& evaluator)
  {
    typename FEFunction::Space::Grid::template Codim<0>::Entity cell = intersection.inside(); //it is important to assign the temporary object to a local variable since the evaluator stores the address of the cell (as pointer) and dereferences it later
    evaluator.moveTo(cell);
    typename FEFunction::ValueType integral(0.0);
    auto const& refElem = evaluator.shapeFunctions().referenceElement();
    auto const& qr = Dune::QuadratureRules<typename FEFunction::Space::ctype,FEFunction::Space::dim-1>::rule(intersection.type(),function.order(evaluator));
    int nQuadPos = qr.size();
    for (int i=0; i<nQuadPos; ++i)
    {
      auto localPositionInCell = refElem.template geometry<1>(intersection.indexInInside()).global(qr[i].position());
      evaluator.evaluateAt(localPositionInCell,qr,i,intersection.indexInInside());
      integral += qr[i].weight()*intersection.geometry().integrationElement(qr[i].position())*function.value(evaluator);
    }
    return integral;
  }
  template <class FEFunction>
  typename FEFunction::ValueType integrateOverBoundary(FEFunction const& function)
  {
    ShapeFunctionCache<typename FEFunction::Space::Grid, typename FEFunction::Space::Scalar> sfCache;
    typename FEFunction::Space::Evaluator evaluator(function.space(),&sfCache);
    typename FEFunction::Space::GridView const& gridView = function.space().gridView();
    typename FEFunction::ValueType integral(0.0);
    forEachBoundaryFace(gridView,[&](typename FEFunction::Space::GridView::Intersection const& intersection)
    {
      integral += integrateOverIntersection(function,intersection,evaluator);
    });
    return integral;
  }
  template <template <class, class> class DomainMapper, typename Scalar, typename GridView, int m>
  typename FunctionSpaceElement<FEFunctionSpace<BoundaryMapper<DomainMapper, Scalar, GridView>>,m>::ValueType
  integrateOverBoundary(FunctionSpaceElement<FEFunctionSpace<BoundaryMapper<DomainMapper, Scalar, GridView>>,m> const& function)
  {
    using FEFunction = FunctionSpaceElement<FEFunctionSpace<BoundaryMapper<DomainMapper, Scalar, GridView>>,m>;
    ShapeFunctionCache<typename GridView::Grid, typename FEFunction::Space::Scalar> sfCache;
    typename FEFunction::Space::Evaluator evaluator(function.space(),&sfCache);
    GridView const& gridView = function.space().gridView();
    BoundaryMapper<DomainMapper, Scalar, GridView> const& mapper = function.space().mapper();
    typename FEFunction::ValueType integral(0.0);
    forEachBoundaryFace(gridView,[&](typename GridView::Intersection const& intersection) {
      if(mapper.inDomain(intersection)) integral += integrateOverIntersection(function,intersection,evaluator);
    });
    return integral;
  }
 
 
  // ---------------------------------------------------------------------------------------------
 
  template <typename Space>
  class Integrator
  {
    typedef typename Space::Grid Grid;
    typedef typename Space::GridView GridView;
 
    typedef typename Grid::template Codim<0>::Entity Cell;
    typedef Dune::QuadraturePoint<typename Grid::ctype,Grid::dimension> QuadPoint;
    typedef typename Space::template Element<1>::type::ValueType ValueType;
 
  public:
    Integrator(): sum(0.0) {}
 
    template <class CellPointer>
    void operator()(CellPointer const& c, QuadPoint const& q, ValueType v)
    {
      v *= q.weight()*c->geometry().integrationElement(q.position());
      sum += v;
    }
 
    ValueType getIntegral() { return sum; }
  private:
    ValueType sum;
  };
 
  template <typename Function,
           class = typename Function::Space>
  typename Function::ValueType integral(Function const& f)
  {
    typename Function::ValueType sum(0.0);
    auto integrator = [&](auto const& cell, auto const& qp, auto const& val)
    {
      sum += qp.weight()*cell.geometry().integrationElement(qp.position())*val;
    };
    forEachQuadPoint(f,integrator);
    return sum;
  }
 
  template <class WeakFunctionView, class Space>
  auto integral(WeakFunctionView const& f, Space const& space)
  {
    Integrator<Space> integrator;
    forEachQuadPoint(f,integrator,space);
    return integrator.getIntegral();
  }
 
  //---------------------------------------------------------------------
  template <typename Space>
  class [[deprecated("Use integral() functions directly. Wil be removed after 2019-12-31.")]] Integral
  {
  public:
    template<typename Function>
    typename Function::ValueType operator()(Function const& f)
    {
      return integral(f);
    }
 
    template<typename Function>
    typename Function::ValueType operator()(Function const& f, typename Function::Space const& space)
    {
      return integral(f,space);
    }
  };
 
 
  struct IntegralF
  {
    template<typename Function>
    typename Function::Scalar id(Function const& f) const
    {
      Integral<typename Function::Space> integral;
      return integral(f);
    }
 
    template<typename Function>
    typename Function::Scalar operator()(Function const& f) const
    {
      return id(f);
    }
  };
 
 
} // end of namespace Kaskade
#endif //FEM_INTEGRATION_HH