integrateentity.hh 1.57 KB
Newer Older
Peter Bastian's avatar
Peter Bastian committed
1 2 3 4 5
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_INTEGRATE_ENTITY_HH
#define DUNE_INTEGRATE_ENTITY_HH

6
#include <dune/common/exceptions.hh>
7
#include <dune/geometry/quadraturerules.hh>
Peter Bastian's avatar
Peter Bastian committed
8 9

//! compute integral of function over entity with given order
10 11
template<class Entity, class Function>
double integrateEntity (const Entity &entity, const Function &f, int p)
Peter Bastian's avatar
Peter Bastian committed
12 13
{
  // dimension of the entity
14
  const int dim = Entity::dimension;
Peter Bastian's avatar
Peter Bastian committed
15 16

  // type used for coordinates in the grid
17
  typedef typename Entity::Geometry::ctype ctype;
18 19 20

  // get geometry
  const typename Entity::Geometry geometry = entity.geometry();
Peter Bastian's avatar
Peter Bastian committed
21 22

  // get geometry type
23
  const Dune::GeometryType gt = geometry.type();
Peter Bastian's avatar
Peter Bastian committed
24 25

  // get quadrature rule of order p
26 27
  const Dune::QuadratureRule<ctype,dim>&
  rule = Dune::QuadratureRules<ctype,dim>::rule(gt,p);     /*@\label{ieh:qr}@*/
Peter Bastian's avatar
Peter Bastian committed
28 29 30 31 32

  // ensure that rule has at least the requested order
  if (rule.order()<p)
    DUNE_THROW(Dune::Exception,"order not available");

33
  // compute approximate integral
Peter Bastian's avatar
Peter Bastian committed
34
  double result=0;
35
  for (typename Dune::QuadratureRule<ctype,dim>::const_iterator i=rule.begin(); /*@\label{ieh:for}@*/
Peter Bastian's avatar
Peter Bastian committed
36 37
       i!=rule.end(); ++i)
  {
38
    double fval = f(geometry.global(i->position()));       /*@\label{ieh:fval}@*/
39
    double weight = i->weight();                        /*@\label{ieh:weight}@*/
40
    double detjac = geometry.integrationElement(i->position());       /*@\label{ieh:detjac}@*/
41
    result += fval * weight * detjac;                   /*@\label{ieh:result}@*/
Peter Bastian's avatar
Peter Bastian committed
42 43
  }

44
  // return result
Peter Bastian's avatar
Peter Bastian committed
45 46
  return result;
}
47

Peter Bastian's avatar
Peter Bastian committed
48
#endif