From e855e7455f78b9ea27de460bf6a98623e549f5a8 Mon Sep 17 00:00:00 2001
From: Felix Albrecht <felix.albrecht@uni-muenster.de>
Date: Wed, 14 Sep 2011 11:46:29 +0200
Subject: [PATCH] began discontinuous galerkin example
---
CMakeLists.txt | 4 +
.../assembler/local/codim1/matrix.hh | 157 +++++++++++++
.../assembler/system/unconstrained.hh | 110 +++++++++
.../discreteoperator/local/codim1/integral.hh | 153 +++++++++++++
.../evaluation/binary/jumpmeanpenalty.hh | 113 ++++++++++
.../elliptic_discontinuous_galerkin/main.cc | 210 ++++++++++++++++++
6 files changed, 747 insertions(+)
create mode 100644 dune/functionals/assembler/local/codim1/matrix.hh
create mode 100644 dune/functionals/assembler/system/unconstrained.hh
create mode 100644 dune/functionals/discreteoperator/local/codim1/integral.hh
create mode 100644 dune/functionals/evaluation/binary/jumpmeanpenalty.hh
create mode 100644 examples/elliptic_discontinuous_galerkin/main.cc
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7cd0adfb9..86aee0515 100644
--- a/CMakeLists.txt
+++ b/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} )
diff --git a/dune/functionals/assembler/local/codim1/matrix.hh b/dune/functionals/assembler/local/codim1/matrix.hh
new file mode 100644
index 000000000..7ce4c2267
--- /dev/null
+++ b/dune/functionals/assembler/local/codim1/matrix.hh
@@ -0,0 +1,157 @@
+#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
diff --git a/dune/functionals/assembler/system/unconstrained.hh b/dune/functionals/assembler/system/unconstrained.hh
new file mode 100644
index 000000000..f6db6ab19
--- /dev/null
+++ b/dune/functionals/assembler/system/unconstrained.hh
@@ -0,0 +1,110 @@
+#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
diff --git a/dune/functionals/discreteoperator/local/codim1/integral.hh b/dune/functionals/discreteoperator/local/codim1/integral.hh
new file mode 100644
index 000000000..95189dc82
--- /dev/null
+++ b/dune/functionals/discreteoperator/local/codim1/integral.hh
@@ -0,0 +1,153 @@
+/**
+ \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
diff --git a/dune/functionals/evaluation/binary/jumpmeanpenalty.hh b/dune/functionals/evaluation/binary/jumpmeanpenalty.hh
new file mode 100644
index 000000000..b76d84abd
--- /dev/null
+++ b/dune/functionals/evaluation/binary/jumpmeanpenalty.hh
@@ -0,0 +1,113 @@
+#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
diff --git a/examples/elliptic_discontinuous_galerkin/main.cc b/examples/elliptic_discontinuous_galerkin/main.cc
new file mode 100644
index 000000000..5d0b80ca6
--- /dev/null
+++ b/examples/elliptic_discontinuous_galerkin/main.cc
@@ -0,0 +1,210 @@
+/**
+ \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;
+ }
+}
--
GitLab