began discontinuous galerkin example

CMakeLists.txt

@@ -127,5 +127,9 @@ execute_process(COMMAND ln -s ${CMAKE_CURRENT_SOURCE_DIR}/examples/macrogrids .)
 ADD_EXECUTABLE( examples/elliptic_finite_element_example "examples/elliptic_finite_element/main.cc" ${common} ${grid} ${istl} ${fem} ${femhowto} ${femtools} ${functionals} )
 TARGET_LINK_LIBRARIES( examples/elliptic_finite_element_example ${COMMON_LIBS}  )
 
+ADD_EXECUTABLE( examples/elliptic_discontinuous_galerkin "examples/elliptic_discontinuous_galerkin/main.cc" ${common} ${grid} ${grid} ${fem} ${femhowto} ${femtools} ${functionals} )
+TARGET_LINK_LIBRARIES( examples/elliptic_discontinuous_galerkin ${COMMON_LIBS} )
+
+
 
 #HEADERCHECK( ${header} )

dune/functionals/assembler/local/codim1/matrix.hh

+#ifndef DUNE_FUNCTIONALS_ASSEMLBER_LOCAL_CODIM1_MATRIX_HH
+#define DUNE_FUNCTIONALS_ASSEMLBER_LOCAL_CODIM1_MATRIX_HH
+
+// local includes
+#include "vector.hh"
+
+namespace Dune {
+
+namespace Functionals {
+
+namespace Assembler {
+
+namespace Local {
+
+namespace Codim1 {
+
+template <class InnerLocalOperatorImp, class DirichletLocalOperatorImp, class NeumannLocalOperatorImp>
+class Matrix
+{
+public:
+  typedef InnerLocalOperatorImp InnerLocalOperatorType;
+
+  typedef DirichletLocalOperatorImp DirichletLocalOperatorType;
+
+  typedef NeumannLocalOperatorImp NeumannLocalOperatorType;
+
+  typedef Matrix<InnerLocalOperatorType, DirichletLocalOperatorType, NeumannLocalOperatorType> ThisType;
+
+  typedef typename InnerLocalOperatorType::RangeFieldType RangeFieldType;
+
+  //  template< class InducingDiscreteFunctionType >
+  //  class LocalVectorAssembler
+  //  {
+  //  private:
+  //    typedef typename LocalOperatorType::template LocalFunctional< InducingDiscreteFunctionType >::Type
+  //      InducingFunctionalType;
+
+  //  public:
+  //    typedef Dune::Functionals::Assembler::Local::Codim0::Vector< InducingFunctionalType >
+  //      Type;
+  //  };
+
+  //! constructor
+  Matrix(const InnerLocalOperatorType& innerLocalOperator, const DirichletLocalOperatorType& dirichletLocalOperator,
+         const NeumannLocalOperatorType& neumannLocalOperator)
+    : innerLocalOperator_(innerLocalOperator)
+    , dirichletLocalOperator_(dirichletLocalOperator)
+    , neumannLocalOperator_(neumannLocalOperator)
+  {
+  }
+
+private:
+  //! copy constructor
+  Matrix(const ThisType& other)
+    : localOperator_(other.localOperator())
+  {
+  }
+
+public:
+  const LocalOperatorType& localOperator() const
+  {
+    return localOperator_;
+  }
+
+  //  template< class InducingDiscreteFunctionType >
+  //  typename LocalVectorAssembler< InducingDiscreteFunctionType >::Type
+  //    localVectorAssembler( const InducingDiscreteFunctionType& inducingDiscreteFunction ) const
+  //  {
+  //    typedef typename LocalVectorAssembler< InducingDiscreteFunctionType >::Type
+  //      LocalVectorAssemblerType;
+
+  //    return LocalVectorAssemblerType( localOperator_.localFunctional( inducingDiscreteFunction ) );
+  //  }
+
+  template <class AnsatzSpaceType, class TestSpaceType, class EntityType, class MatrixType, class LocalMatrixType>
+  void assembleLocal(const AnsatzSpaceType& ansatzSpace, const TestSpaceType& testSpace, const EntityType& entity,
+                     MatrixType& matrix, LocalMatrixType& tmpLocalMatrix) const
+  {
+    typedef typename AnsatzSpaceType::GridPartType GridPartType;
+
+    typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
+
+    typedef typename IntersectionIteratorType::Intersection IntersectionType;
+
+    typedef typename IntersectionType::EntityPointer EntityPointerType;
+
+    const GridPartType& gridPart = ansatzSpace.gridPart();
+
+    const IntersectionIteratorType lastIntersection = gridPart.iend(entity);
+
+    for (IntersectionIteratorType intIt = entity.ibegin(); intIt != lastIntersection; ++intIt) {
+      const IntersectionType& intersection = *intIt;
+
+      // if inner intersection
+      if (intersection.neighbor() && !intersection.boundary()) {
+        const EntityPointerType neighbourPtr = intersection.outside();
+        const EntityType& neighbour          = *neighbourPtr;
+
+        innerLocalOperator_.applyLocal(intersection,
+                                       ansatzSpace.baseFunctionSet().local(entity),
+                                       testSpace.baseFunctionSet().local(neighbour),
+                                       tmpLocalMatrix);
+
+        // write local matrix to global
+        addToMatrix(ansatzSpace, testSpace, entity, neighbour, tmpLocalMatrix, matrix);
+
+
+      } else if (!intersection.neighbor() && intersection.boundary()) // else boundary
+      {
+        const unsigned int boundaryId = intersection.boundaryId();
+
+        // if dirichlet
+        if (boundaryId == 2) {
+
+        } else if (boundaryId == 3) // else neumann
+        {
+        }
+
+      } // end if boundary intersection
+    }
+  }
+
+private:
+  //! assignment operator
+  ThisType& operator=(const ThisType&);
+
+  template <class AnsatzSpaceType, class TestSpaceType, class EntityType, class LocalMatrixType, class MatrixType>
+  void addToMatrix(const AnsatzSpaceType& ansatzSpace, const TestSpaceType& testSpace, const EntityType& entity,
+                   const EntityType& neighbour, const LocalMatrixType& localMatrix, MatrixType& matrix) const
+  {
+    for (unsigned int i = 0; i < ansatzSpace.baseFunctionSet().local(entity).size(); ++i) {
+      for (unsigned int j = 0; j < testSpace.baseFunctionSet().local(entity).size(); ++j) {
+        const unsigned int globalI = ansatzSpace.map().toGlobal(entity, i);
+        const unsigned int globalJ = testSpace.map().toGlobal(neighbour, j);
+
+        matrix[globalI][globalJ] += localMatrix[i][j];
+      }
+    }
+  } // end method addToMatrix
+
+  const InnerLocalOperatorType& innerLocalOperator_;
+  const DirichletLocalOperatorType& dirichletLocalOperator_;
+  const NeumannLocalOperatorType& neumannLocalOperator_;
+
+}; // end class Matrix
+
+} // end namespace Codim1
+
+} // end namespace Local
+
+} // end namespace Assembler
+
+} // end namespace Functionals
+
+} // end namespace Dune
+
+#endif // DUNE_FUNCTIONALS_ASSEMLBER_LOCAL_CODIM1_MATRIX_HH

dune/functionals/assembler/system/unconstrained.hh

+#ifndef DUNE_FUNCTIONALS_ASSEMBLER_SYSTEM_UNCONSTRAINED_HH
+#define DUNE_FUNCTIONALS_ASSEMBLER_SYSTEM_UNCONSTRAINED_HH
+
+// dune-common includes
+#include <dune/common/dynmatrix.hh>
+#include <dune/common/dynvector.hh>
+
+namespace Dune {
+
+namespace Functionals {
+
+namespace Assembler {
+
+namespace System {
+
+template <class AnsatzFunctionSpaceImp, class TestFunctionSpaceImp = AnsatzFunctionSpaceImp>
+class Unconstrained
+{
+public:
+  typedef AnsatzFunctionSpaceImp AnsatzFunctionSpaceType;
+
+  typedef TestFunctionSpaceImp TestFunctionSpaceType;
+
+  typedef Unconstrained<AnsatzFunctionSpaceImp, TestFunctionSpaceImp> ThisType;
+
+  //! constructor
+  Unconstrained(const AnsatzFunctionSpaceType& ansatzSpace, const TestFunctionSpaceType& testSpace)
+    : ansatzSpace_(ansatzSpace)
+    , testSpace_(testSpace)
+  {
+  }
+
+  //! constructor
+  Unconstrained(const AnsatzFunctionSpaceType& ansatzSpace)
+    : ansatzSpace_(ansatzSpace)
+    , testSpace_(ansatzSpace)
+  {
+  }
+
+private:
+  //! copy constructor
+  Unconstrained(const ThisType& other)
+    : ansatzSpace_(other.ansatzSpace())
+    , testSpace_(other.testSpace())
+  {
+  }
+
+public:
+  const AnsatzFunctionSpaceType& ansatzSpace()
+  {
+    return ansatzSpace_;
+  }
+
+  const TestFunctionSpaceType& testSpace()
+  {
+    return testSpace_;
+  }
+
+  template <class LocalMatrixAssemblerType, class MatrixType, class LocalVectorAssemblerType, class VectorType>
+  void assembleSystem(const LocalMatrixAssemblerType& localMatrixAssembler, MatrixType& systemMatrix,
+                      const LocalVectorAssemblerType& localVectorAssembler, VectorType& systemVector) const
+  {
+    // some types
+    typedef typename AnsatzFunctionSpaceType::GridPartType GridPartType;
+
+    typedef typename AnsatzFunctionSpaceType::IteratorType EntityIteratorType;
+
+    typedef typename AnsatzFunctionSpaceType::EntityType EntityType;
+
+    typedef typename AnsatzFunctionSpaceType::RangeFieldType RangeFieldType;
+
+    typedef Dune::DynamicMatrix<RangeFieldType> LocalMatrixType;
+
+    typedef Dune::DynamicVector<RangeFieldType> LocalVectorType;
+
+    // common storage for all entities
+    LocalMatrixType tmpLocalMatrix(
+        ansatzSpace_.map().maxLocalSize(), testSpace_.map().maxLocalSize(), RangeFieldType(0.0));
+    LocalVectorType tmpLocalVector(testSpace_.map().maxLocalSize(), RangeFieldType(0.0));
+
+    // do gridwalk to assemble
+    const EntityIteratorType lastEntity = ansatzSpace_.end();
+    for (EntityIteratorType entityIterator = ansatzSpace_.begin(); entityIterator != lastEntity; ++entityIterator) {
+      const EntityType& entity = *entityIterator;
+
+      localMatrixAssembler.assembleLocal(ansatzSpace_, testSpace_, entity, systemMatrix, tmpLocalMatrix);
+      localVectorAssembler.assembleLocal(testSpace_, entity, systemVector, tmpLocalVector);
+
+    } // done gridwalk to assemble
+
+  } // end method assembleSystem
+
+private:
+  //! assignment operator
+  ThisType& operator=(const ThisType&);
+
+  const AnsatzFunctionSpaceType& ansatzSpace_;
+  const TestFunctionSpaceType& testSpace_;
+
+}; // end class Unconstrained
+
+} // end namespace System
+
+} // end namespace Assembler
+
+} // end namespace Functionals
+
+} // end namespace Dune
+
+#endif // DUNE_FUNCTIONALS_ASSEMBLER_SYSTEM_UNCONSTRAINED_HH

dune/functionals/discreteoperator/local/codim1/integral.hh

+/**
+  \file   integral.hh
+  **/
+
+#ifndef DUNE_FUNCTIONALS_DISCRETEOPERATOR_LOCAL_CODIM1_INTEGRAL_HH
+#define DUNE_FUNCTIONALS_DISCRETEOPERATOR_LOCAL_CODIM1_INTEGRAL_HH
+
+// dune fem includes
+#include <dune/fem/quadrature/cachingquadrature.hh>
+
+// dune-functionals includes
+//#include <dune/functionals/common/localmatrix.hh>
+//#include <dune/functionals/common/localvector.hh>
+//#include <dune/functionals/discretefunctional/local/codim0/integral.hh>
+
+namespace Dune {
+
+namespace Functionals {
+
+namespace DiscreteOperator {
+
+namespace Local {
+
+namespace Codim1 {
+
+template <class LocalEvaluationImp>
+class Integral
+{
+public:
+  typedef LocalEvaluationImp LocalEvaluationType;
+
+  typedef Integral<LocalEvaluationType> ThisType;
+
+  typedef typename LocalEvaluationType::FunctionSpaceType FunctionSpaceType;
+
+  typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
+
+  typedef typename FunctionSpaceType::DomainType DomainType;
+
+  //  template< class InducingDiscreteFunctionType >
+  //  class LocalFunctional
+  //  {
+  //  public:
+  //    typedef Dune::Functionals::DiscreteFunctional::Local::Codim0::IntegralInduced<  ThisType,
+  //                                                                                    InducingDiscreteFunctionType >
+  //      Type;
+  //  };
+
+  Integral(const LocalEvaluationType localEvaluation)
+    : localEvaluation_(localEvaluation)
+  {
+  }
+
+  //! copy constructor
+  Integral(const ThisType& other)
+    : localEvaluation_(other.localEvaluation())
+  {
+  }
+
+  const LocalEvaluationType& localEvaluation() const
+  {
+    return localEvaluation_;
+  }
+
+  //  template< class InducingDiscreteFunctionType >
+  //  typename LocalFunctional< InducingDiscreteFunctionType >::Type
+  //    localFunctional( const InducingDiscreteFunctionType& inducingDiscreteFunction ) const
+  //  {
+  //    typedef Dune::Functionals::DiscreteFunctional::Local::Codim0::IntegralInduced<  ThisType,
+  //                                                                                    InducingDiscreteFunctionType >
+  //      LocalFunctionalType;
+
+  //    return LocalFunctionalType( *this, inducingDiscreteFunction );
+  //  } // end method localFunctional
+
+  template <class IntersectionType, class LocalAnsatzBaseFunctionSetType, /*Entity*/
+            class LocalTestBaseFunctionSetType, /*Neighour*/
+            class LocalMatrixType>
+  void applyLocal(const IntersectionType& intersection,
+                  const LocalAnsatzBaseFunctionSetType& localAnsatzBaseFunctionSet,
+                  const LocalTestBaseFunctionSetType& localTestBaseFunctionSet, LocalMatrixType& localMatrix) const
+  {
+    // clear target matrix
+    for (unsigned int i = 0; i < localMatrix.rows(); ++i) {
+      for (unsigned int j = 0; j < localMatrix.cols(); ++j) {
+        localMatrix[i][j] = 0.0;
+      }
+    }
+
+    // some types
+    typedef typename LocalAnsatzBaseFunctionSetType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
+
+    typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
+
+    typedef Dune::CachingQuadrature<GridPartType, 1> FaceQuadratureType;
+
+    // some stuff
+    const GridPartType& gridPart = localAnsatzBaseFunctionSet.space().gridPart();
+    const unsigned int rows      = localAnsatzBaseFunctionSet.size();
+    const unsigned int cols = localTestBaseFunctionSet.size();
+    const unsigned int quadratureOrder =
+        localEvaluation_.order() + localAnsatzBaseFunctionSet.order() + localTestBaseFunctionSet.order();
+    const FaceQuadratureType faceQuadrature(gridPart, intersection, quadratureOrder, FaceQuadratureType::INSIDE);
+    const unsigned int numberOfQuadraturePoints = faceQuadrature.nop();
+
+    // some tmp storage
+    LocalMatrixType tmpMatrix(rows, cols);
+
+    // do loop over all quadrature points
+    for (unsigned int q = 0; q < numberOfQuadraturePoints; ++q) {
+      // local coordinates
+      const DomainType x               = faceQuadrature.point(q);
+      const DomainType xInEntity       = intersection.geometryInInside().global(x);
+      const DomainType xInNeighbour    = intersection.geometryInOutside().global(x);
+      const DomainType unitOuterNormal = intersection.unitOuterNormal();
+
+      // integration factors
+      const double integrationFactor = intersection.geometry().integrationElement(x);
+      const double quadratureWeight  = faceQuadrature.weight(q);
+
+      // evaluate the local operation
+      localEvaluation_.evaluate(
+          localAnsatzBaseFunctionSet, localTestBaseFunctionSet, xInEntity, xInNeighbour, unitOuterNormal, tmpMatrix);
+
+      // compute integral
+      for (unsigned int i = 0; i < rows; ++i) {
+        for (unsigned int j = 0; j < cols; ++j) {
+          localMatrix[i][j] += tmpMatrix[i][j] * integrationFactor * quadratureWeight;
+        }
+      }
+    } // done loop over all quadrature points
+
+  } // end method applyLocal
+
+private:
+  //! assignment operator
+  ThisType& operator=(const ThisType&);
+
+  const LocalEvaluationType localEvaluation_;
+
+}; // end class Codim0Integration
+
+} // end namespace Codim1
+
+} // end namespace Local
+
+} // end namespace DiscreteOperator
+
+} // end namespace Functionals
+
+} // end namespace Dune
+
+#endif // end DUNE_FUNCTIONALS_DISCRETEOPERATOR_LOCAL_CODIM1_INTEGRAL_HH

dune/functionals/evaluation/binary/jumpmeanpenalty.hh

+#ifndef DUNE_FUNCTIONALS_EVALUATION_BINARY_JUMPMEANPENALTY_HH
+#define DUNE_FUNCTIONALS_EVALUATION_BINARY_JUMPMEANPENALTY_HH
+
+namespace Dune {
+
+namespace Functionals {
+
+namespace Evaluation {
+
+namespace Binary {
+
+/**
+  \brief  This represents the operation \f$a\nabla u \nabla v\f$.
+
+          \f$a\f$ is a given scalar function (in this case 1) and \f$u\f$ and \f$u\f$ may be local functions, i.e.
+          ansatz- and testfunctions.
+  \tparam FunctionSpaceImp
+          Type of the function space, where \f$f\f$, \f$u\f$ and \f$v\f$ live in.
+  **/
+template <class FunctionSpaceImp>
+class JumpMeanPenalty
+{
+public:
+  typedef FunctionSpaceImp FunctionSpaceType;
+
+  typedef JumpMeanPenalty<FunctionSpaceType> ThisType;
+
+  typedef typename FunctionSpaceType::DomainType DomainType;
+
+  typedef typename FunctionSpaceType::RangeType RangeType;
+
+  typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
+
+  typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
+
+  typedef Dune::FemTools::Function::Runtime<FunctionSpaceType> InducingFunctionType;
+
+  //! constructor, takes the inducing functions expression as a runtime parameter
+  JumpMeanPenalty(const std::string expression = "[1.0;1.0;1.0]", const int order = 1)
+    : inducingFunction_(expression)
+    , order_(std::max(0, order))
+  {
+  }
+
+  //! copy constructor
+  JumpMeanPenalty(const ThisType& other)
+    : inducingFunction_(other.inducingFunction())
+    , order_(other.order())
+  {
+  }
+
+  //! returns the inducing function
+  InducingFunctionType inducingFunction() const
+  {
+    return inducingFunction_;
+  }
+
+  unsigned int order() const
+  {
+    return order_;
+  }
+
+  template <class LocalAnsatzBaseFunctionSetType, class LocalTestBaseFunctionSetType, class LocalMatrixType>
+  void evaluate(const LocalAnsatzBaseFunctionSetType& localAnsatzBaseFunctionSet,
+                const LocalTestBaseFunctionSetType& localTestBaseFunctionSet, const DomainType& localPointEntity,
+                const DomainType& localPointNeighbour, const DomainType& unitOuterNormal, LocalMatrixType& ret) const
+  {
+    // get global point
+    const DomainType globalPoint = localAnsatzBaseFunctionSet.entity().geometry().global(localPointEntity);
+
+    // evaluate first gradient
+    const unsigned int rows = localAnsatzBaseFunctionSet.size();
+    std::vector<JacobianRangeType> gradientLocalAnsatzBaseFunctionSet(rows, JacobianRangeType(0.0));
+    localAnsatzBaseFunctionSet.jacobian(localPointEntity, gradientLocalAnsatzBaseFunctionSet);
+
+    // evaluate second gradient
+    const unsigned int cols = localTestBaseFunctionSet.size();
+    std::vector<JacobianRangeType> gradientLocalTestBaseFunctionSet(cols, JacobianRangeType(0.0));
+    localTestBaseFunctionSet.jacobian(localPointNeighbour, gradientLocalTestBaseFunctionSet);
+
+    // evaluate inducing function
+    RangeType functionValue(0.0);
+    inducingFunction_.evaluate(globalPoint, functionValue);
+
+    // do loop over all ansatz and test basefunctions
+    assert(ret.rows() == rows);
+    assert(ret.cols() == cols);
+    for (unsigned int i = 0; i < rows; ++i) {
+      for (unsigned int j = 0; j < cols; ++j) {
+        const RangeFieldType gradientProduct =
+            gradientLocalAnsatzBaseFunctionSet[i][0] * gradientLocalTestBaseFunctionSet[j][0];
+        ret[i][j] = functionValue * gradientProduct;
+      }
+    }
+  }
+
+private:
+  //! assignment operator
+  ThisType& operator=(const ThisType&);
+
+  const InducingFunctionType inducingFunction_;
+  unsigned int order_;
+}; // end class JumpMeanPenalty
+
+} // end namespace Binary
+
+} // end namespace Evaluation
+
+} // end namespace Functionals
+
+} // end namespace Dune
+
+#endif // DUNE_FUNCTIONALS_EVALUATION_BINARY_JUMPMEANPENALTY_HH

examples/elliptic_discontinuous_galerkin/main.cc

+/**
+  \file   main.cc
+  \brief  Main file for the elliptic discontinuous galerkin example.
+  **/
+
+// disable warnings about problems in dune headers
+#include <dune/fem-tools/header/disablewarnings.hh>
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <iostream>
+#include <vector>
+
+// dune-common includes
+#include <dune/common/exceptions.hh>
+//#include <dune/common/fmatrix.hh>
+
+// dune-grid includes
+#include <dune/grid/utility/gridtype.hh>
+
+// dune-istl includes
+#include <dune/istl/solvers.hh>
+
+// dune-fem includes
+//#include <dune/fem/function/adaptivefunction.hh>
+//#include <dune/fem/function/blockvectorfunction.hh>
+#include <dune/fem/gridpart/gridpart.hh>
+//#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
+#include <dune/fem/misc/mpimanager.hh>
+
+// reenable warnings
+#include <dune/fem-tools/header/enablewarnings.hh>
+
+// dune-functionals includes
+#include <dune/functionals/discretefunctionspace/discontinuous/orthonormal.hh>
+//#include <dune/functionals/container/factory.hh>
+//#include <dune/functionals/assembler/local/codim0/matrix.hh>
+//#include <dune/functionals/assembler/local/codim0/vector.hh>
+//#include <dune/functionals/assembler/system/unconstrained.hh>
+//#include <dune/functionals/discretefunction/continuous.hh>
+//#include <dune/functionals/discretefunction/femadapter.hh>
+//#include <dune/functionals/discreteoperator/local/codim0/integral.hh>
+//#include <dune/functionals/discretefunctional/local/codim0/integral.hh>
+//#include <dune/functionals/evaluation/unary/scale.hh>
+//#include <dune/functionals/evaluation/binary/elliptic.hh>
+
+// dune-fem-tools includes
+#include <dune/fem-tools/common/string.hh>
+#include <dune/fem-tools/common/printing.hh>
+//#include <dune/fem-tools/function/runtimefunction.hh>
+#include <dune/fem-tools/discretefunction.hh>
+
+using namespace Dune::Functionals;
+
+#ifndef POLORDER
+const int polOrder = 1;
+#else
+const int polOrder = POLORDER;
+#endif
+
+
+int main(int argc, char** argv)
+{
+  try {
+
+    // MPI manager
+    Dune::MPIManager::initialize(argc, argv);
+
+
+    // grid
+    typedef Dune::GridSelector::GridType GridType;
+
+    typedef Dune::LeafGridPart<GridType> GridPartType;
+
+    const std::string dgfFilename = "../macrogrids/unitcube" + Dune::FemTools::String::toString(GRIDDIM) + ".dgf";
+
+    Dune::GridPtr<GridType> gridPtr(dgfFilename);
+
+    GridPartType gridPart(*gridPtr);
+
+    static const unsigned int dimDomain = GridType::dimension;
+
+    static const unsigned int dimRange = 1;
+
+    // function space
+    typedef Dune::FunctionSpace<double, double, dimDomain, dimRange> FunctionSpaceType;
+
+    typedef /*typename*/ FunctionSpaceType::RangeFieldType RangeFieldType;
+
+
+    // discrete function space
+    typedef DiscreteFunctionSpace::Discontinuous::Orthonormal<FunctionSpaceType, GridPartType, polOrder> DGSpaceType;
+
+    const DGSpaceType dgSpace(gridPart);
+
+
+    //    // local evaluation
+    //    typedef Dune::Functionals::Evaluation::Unary::Scale< FunctionSpaceType >
+    //      ProductEvaluationType;
+
+    //    ProductEvaluationType productEvaluation( "[1.0;1.0;1.0]", 0 );
+
+    //    typedef Dune::Functionals::Evaluation::Binary::Elliptic< FunctionSpaceType >
+    //      EllipticEvaluationType;
+
+    //    EllipticEvaluationType ellipticEvaluation( "[1.0;1.0;1.0]", 0 );
+
+
+    //    // operator and functional
+    //    typedef DiscreteOperator::Local::Codim0::Integral< EllipticEvaluationType >
+    //      LocalEllipticOperatorType;
+
+    //    const LocalEllipticOperatorType localEllipticOperator( ellipticEvaluation );
+
+    //    typedef DiscreteFunctional::Local::Codim0::Integral< ProductEvaluationType >
+    //      LocalL2FunctionalType;
+
+    //    const LocalL2FunctionalType localL2Functional( productEvaluation );
+
+
+    //    // matrix, rhs and solution storage
+    //    typedef Container::Matrix::Defaults< RangeFieldType, dimRange, dimRange >::BCRSMatrix
+    //      MatrixFactory;
+
+    //    typedef /*typename*/ MatrixFactory::AutoPtrType
+    //      MatrixPtrType;
+
+    //    MatrixPtrType A = MatrixFactory::create( discreteH1, discreteH1 );
+
+    //    typedef Container::Vector::Defaults< RangeFieldType, dimRange >::BlockVector
+    //      VectorFactory;
+
+    //    typedef /*typename*/ VectorFactory::AutoPtrType
+    //      VectorPtrType;
+
+    //    VectorPtrType F = VectorFactory::create( discreteH1 );
+    //    *F = 0.0;
+
+    //    VectorPtrType u = VectorFactory::create( discreteH1 );
+    //    *u = 0.0;
+
+
+    //    // assembler
+    //    typedef Assembler::Local::Codim0::Matrix< LocalEllipticOperatorType >
+    //      LocalMatrixAssemblerType;
+
+    //    const LocalMatrixAssemblerType localMatrixAssembler( localEllipticOperator );
+
+    //    typedef Assembler::Local::Codim0::Vector< LocalL2FunctionalType >
+    //      LocalVectorAssemblerType;
+
+    //    const LocalVectorAssemblerType localVectorAssembler( localL2Functional );
+
+    //    typedef Assembler::System::Unconstrained< DiscreteH1Type, DiscreteH1Type >
+    //      SystemAssemblerType;
+
+    //    const SystemAssemblerType systemAssembler( discreteH1, discreteH1 );
+
+    //    systemAssembler.assembleSystem( localMatrixAssembler, *A,
+    //                                    localVectorAssembler, *F );
+
+
+    //    // preconditioner and solver
+    //    typedef /*typename*/ MatrixFactory::ContainerType
+    //      MatrixContainerType;
+
+    //    typedef /*typename*/ VectorFactory::ContainerType
+    //      VectorContainerType;
+
+    //    typedef Dune::MatrixAdapter< MatrixContainerType, VectorContainerType, VectorContainerType >
+    //      MatrixAdapterType;
+
+    //    MatrixAdapterType matrix( *A );
+
+    //    typedef Dune::SeqILU0< MatrixContainerType, VectorContainerType, VectorContainerType, 1 >
+    //      PreconditionerType;
+
+    //    PreconditionerType preconditioner( *A, 1.0 );
+
+    //    typedef Dune::CGSolver< VectorContainerType >
+    //      SolverType;
+
+    //    SolverType solver( matrix, preconditioner, 1e-4, 100, 2 );
+
+    //    Dune::InverseOperatorResult result;
+
+    //    // u_0 = A^(-1) F
+    //    solver.apply( *u, *F, result );
+
+
+    //    // postprocessing
+    //    typedef /*typename*/ Dune::Functionals::DiscreteFunction::Continuous::BlockVector< DiscreteH1Type >
+    //      DiscreteFunctionType;
+
+    //    DiscreteFunctionType solution( discreteH1, *u, "solution" );
+    //    Dune::FemTools::DiscreteFunction::IO::writeToVTK( solution, "solution" );
+
+
+    // done
+    return 0;
+  } catch (Dune::Exception& e) {
+    std::cerr << "Dune reported error: " << e << std::endl;
+  } catch (std::exception& e) {
+    std::cerr << e.what() << std::endl;
+  } catch (...) {
+    std::cerr << "Unknown exception thrown!" << std::endl;
+  }
+}