diff --git a/configure.ac b/configure.ac
index eb5295eb9ac77aa38fc4de3d0d5207a6a3eda98c..131c10ff78eb0ec13f7706b0cf7683a2332621ad 100644
--- a/configure.ac
+++ b/configure.ac
@@ -27,6 +27,7 @@ AC_CONFIG_FILES([
dune/detailed-discretizations/assembler/Makefile
dune/detailed-discretizations/assembler/local/Makefile
dune/detailed-discretizations/assembler/local/codim0/Makefile
+ dune/detailed-discretizations/assembler/local/codim1/Makefile
dune/detailed-discretizations/assembler/system/Makefile
dune/detailed-discretizations/basefunctionset/Makefile
dune/detailed-discretizations/basefunctionset/continuous/Makefile
@@ -38,15 +39,18 @@ AC_CONFIG_FILES([
dune/detailed-discretizations/discretefunctional/Makefile
dune/detailed-discretizations/discretefunctional/local/Makefile
dune/detailed-discretizations/discretefunctional/local/codim0/Makefile
+ dune/detailed-discretizations/discretefunctional/local/codim1/Makefile
dune/detailed-discretizations/discretefunctionspace/Makefile
dune/detailed-discretizations/discretefunctionspace/continuous/Makefile
dune/detailed-discretizations/discretefunctionspace/subspace/Makefile
dune/detailed-discretizations/discreteoperator/Makefile
dune/detailed-discretizations/discreteoperator/local/Makefile
dune/detailed-discretizations/discreteoperator/local/codim0/Makefile
+ dune/detailed-discretizations/discreteoperator/local/codim1/Makefile
dune/detailed-discretizations/evaluation/Makefile
dune/detailed-discretizations/evaluation/local/Makefile
dune/detailed-discretizations/evaluation/local/binary/Makefile
+ dune/detailed-discretizations/evaluation/local/quaternary/Makefile
dune/detailed-discretizations/evaluation/local/unary/Makefile
dune/detailed-discretizations/mapper/Makefile
dune/detailed-discretizations/mapper/continuous/Makefile
diff --git a/dune/detailed-discretizations/assembler/local/Makefile.am b/dune/detailed-discretizations/assembler/local/Makefile.am
index d931f05f8bee4396e467ddedc6e45a67a87fe4c1..827ea7f1590c698b75acb3956807b61eb58794e7 100644
--- a/dune/detailed-discretizations/assembler/local/Makefile.am
+++ b/dune/detailed-discretizations/assembler/local/Makefile.am
@@ -1,3 +1,3 @@
-SUBDIRS = codim0
+SUBDIRS = codim0 codim1
include $(top_srcdir)/am/global-rules
diff --git a/dune/detailed-discretizations/assembler/local/codim1/Makefile.am b/dune/detailed-discretizations/assembler/local/codim1/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..819957a4e6bf67fc8dd3b386a31649ee012a631a
--- /dev/null
+++ b/dune/detailed-discretizations/assembler/local/codim1/Makefile.am
@@ -0,0 +1,4 @@
+assembler_local_codim1_includedir = $(includedir)/dune/detailed-discretizations/assembler/local/codim1
+assembler_local_codim1_include_HEADERS = matrix.hh
+
+include $(top_srcdir)/am/global-rules
diff --git a/dune/detailed-discretizations/assembler/local/codim1/matrix.hh b/dune/detailed-discretizations/assembler/local/codim1/matrix.hh
index 5bfe0ded019375cb9c0af5bf6e1bfd6e8d9ea77b..eeb8202bbf2c2168b2eaf398a90645e7b18c833e 100644
--- a/dune/detailed-discretizations/assembler/local/codim1/matrix.hh
+++ b/dune/detailed-discretizations/assembler/local/codim1/matrix.hh
@@ -17,19 +17,26 @@ namespace Local {
namespace Codim1 {
-template <class InnerLocalOperatorImp, class DirichletLocalOperatorImp, class NeumannLocalOperatorImp>
+/**
+ \todo Add neumann boundary treatment
+ \todo When adding neumann: think of numTmpObjectsRequired()!
+ \todo Add vector assembler
+ \todo Add penalty parameter
+ **/
+template <class LocalInnerOperatorImp, class LocalDirichletOperatorImp /*, class LocalNeumannOperatorImp*/>
class Matrix
{
public:
- typedef InnerLocalOperatorImp InnerLocalOperatorType;
+ typedef LocalInnerOperatorImp LocalInnerOperatorType;
- typedef DirichletLocalOperatorImp DirichletLocalOperatorType;
+ typedef LocalDirichletOperatorImp LocalDirichletOperatorType;
- typedef NeumannLocalOperatorImp NeumannLocalOperatorType;
+ // typedef LocalNeumannOperatorImp
+ // LocalNeumannOperatorType;
- typedef Matrix<InnerLocalOperatorType, DirichletLocalOperatorType, NeumannLocalOperatorType> ThisType;
+ typedef Matrix<LocalInnerOperatorType, LocalDirichletOperatorType /*, LocalNeumannOperatorType*/> ThisType;
- typedef typename InnerLocalOperatorType::RangeFieldType RangeFieldType;
+ typedef typename LocalInnerOperatorType::RangeFieldType RangeFieldType;
// template< class InducingDiscreteFunctionType >
// class LocalVectorAssembler
@@ -44,27 +51,39 @@ public:
// };
//! constructor
- Matrix(const InnerLocalOperatorType innerLocalOperator, const DirichletLocalOperatorType dirichletLocalOperator,
- const NeumannLocalOperatorType neumannLocalOperator)
- : innerLocalOperator_(innerLocalOperator)
- , dirichletLocalOperator_(dirichletLocalOperator)
- , neumannLocalOperator_(neumannLocalOperator)
+ Matrix(const LocalInnerOperatorType localInnerOperator, const LocalDirichletOperatorType localDirichletOperator /*,
+ const LocalNeumannOperatorType localNeumannOperator*/)
+ : localInnerOperator_(localInnerOperator)
+ , localDirichletOperator_(localDirichletOperator) /*,
+ localNeumannOperator_( localNeumannOperator )*/
{
}
private:
//! copy constructor
Matrix(const ThisType& other)
- : localOperator_(other.localOperator())
+ : localInnerOperator_(other.localInnerOperator())
+ , localDirichletOperator_(other.localDirichletOperator()) /*,
+ localNeumannOperator_( other.localNeumannOperator() )*/
{
}
public:
- const LocalOperatorType& localOperator() const
+ const LocalInnerOperatorType& localInnerOperator() const
{
- return localOperator_;
+ return localInnerOperator_;
}
+ const LocalDirichletOperatorType& localDirichletOperator() const
+ {
+ return localDirichletOperator_;
+ }
+
+ // const LocalNeumannOperatorType& localNeumannOperator() const
+ // {
+ // return localNeumannOperator_;
+ // }
+
// template< class InducingDiscreteFunctionType >
// typename LocalVectorAssembler< InducingDiscreteFunctionType >::Type
// localVectorAssembler( const InducingDiscreteFunctionType& inducingDiscreteFunction ) const
@@ -75,23 +94,36 @@ public:
// return LocalVectorAssemblerType( localOperator_.localFunctional( inducingDiscreteFunction ) );
// }
- template <class AnsatzSpaceType, class TestSpaceType, class EntityType, class MatrixType, class LocalMatrixType>
+ /**
+ \todo Add neumann treatment here!
+ **/
+ std::vector<unsigned int> numTmpObjectsRequired() const
+ {
+ std::vector<unsigned int> ret(2, 0);
+ // we require 4 tmp matrix in this local assembler
+ ret[0] = 4;
+ // the operator itself requires that much local matrices
+ ret[1] = std::max(localInnerOperator_.numTmpObjectsRequired(), localDirichletOperator_.numTmpObjectsRequired());
+ return ret;
+ }
+
+ template <class AnsatzSpaceType, class TestSpaceType, class EntityType, class SystemMatrixType, class LocalMatrixType>
void assembleLocal(const AnsatzSpaceType& ansatzSpace, const TestSpaceType& testSpace, const EntityType& entity,
- MatrixType& matrix, std::vector<LocalMatrixType>& tmpLocalMatrices) const
+ SystemMatrixType& systemMatrix,
+ std::vector<std::vector<LocalMatrixType>>& tmpLocalMatricesContainer) const
{
// get the local basefunction sets
- typedef typename AnsatzSpaceType::BaseFunctionSet::Local LocalAnsatzBaseFunctionSetType;
+ typedef typename AnsatzSpaceType::BaseFunctionSetType::LocalBaseFunctionSetType LocalAnsatzBaseFunctionSetType;
- const LocalAnsatzBaseFunctionSetType localAnsatzBaseFunctionSetEntity = ansatzSpace.baseFunctionSet().local(entity);
- const LocalAnsatzBaseFunctionSetType localAnsatzBaseFunctionSetNeighbour =
- ansatzSpace.baseFunctionSet().local(neighbour);
+ const LocalAnsatzBaseFunctionSetType localAnsatzBaseFunctionSetEn = ansatzSpace.baseFunctionSet().local(entity);
- typedef typename TestSpaceType::BaseFunctionSet::Local LocalTesBaseFunctionSetType;
+ typedef typename TestSpaceType::BaseFunctionSetType::LocalBaseFunctionSetType LocalTesBaseFunctionSetType;
- const LocalTesBaseFunctionSetType localTesBaseFunctionSetEntity = testSpace.baseFunctionSet().local(entity);
- const LocalTesBaseFunctionSetType localTesBaseFunctionSetNeighbour = testSpace.baseFunctionSet().local(neighbour);
+ const LocalTesBaseFunctionSetType localTesBaseFunctionSetEn = testSpace.baseFunctionSet().local(entity);
- // check, if we have enough tmp matrices
+ // check tmp local matrices
+ assert(tmpLocalMatricesContainer.size() > 1);
+ std::vector<LocalMatrixType>& tmpLocalMatrices = tmpLocalMatricesContainer[0];
if (tmpLocalMatrices.size() < 4) {
tmpLocalMatrices.resize(
4, LocalMatrixType(ansatzSpace.map().maxLocalSize(), testSpace.map().maxLocalSize(), RangeFieldType(0.0)));
@@ -116,62 +148,86 @@ public:
// if inner intersection
if (intersection.neighbor() && !intersection.boundary()) {
+ // get neighbouring entity
const EntityPointerType neighbourPtr = intersection.outside();
const EntityType& neighbour = *neighbourPtr;
- innerLocalOperator_.applyLocal(localAnsatzBaseFunctionSetEntity,
- localAnsatzBaseFunctionSetNeighbour,
- localTesBaseFunctionSetEntity,
- localTesBaseFunctionSetNeighbour,
- intersection,
- tmpLocalMatrices[0],
- tmpLocalMatrices[1],
- tmpLocalMatrices[2],
- tmpLocalMatrices[3]);
-
- // write local matrix to global
- addToMatrix(ansatzSpace, testSpace, entity, neighbour, tmpLocalMatrix, matrix);
+ // do visit only once
+ if (gridPart.indexSet().index(entity) < gridPart.indexSet().index(neighbour)) {
+ // get neighbouring local basefunction sets
+ const LocalAnsatzBaseFunctionSetType localAnsatzBaseFunctionSetNe =
+ ansatzSpace.baseFunctionSet().local(neighbour);
+ const LocalTesBaseFunctionSetType localTesBaseFunctionSetNe = testSpace.baseFunctionSet().local(neighbour);
+
+ localInnerOperator_.applyLocal(localAnsatzBaseFunctionSetEn,
+ localAnsatzBaseFunctionSetNe,
+ localTesBaseFunctionSetEn,
+ localTesBaseFunctionSetNe,
+ intersection,
+ tmpLocalMatrices[0],
+ tmpLocalMatrices[1],
+ tmpLocalMatrices[2],
+ tmpLocalMatrices[3],
+ tmpLocalMatricesContainer[1]);
+
+ // write local matrix to global (see below)
+ addToMatrix(ansatzSpace, testSpace, entity, entity, tmpLocalMatrices[0], systemMatrix);
+ addToMatrix(ansatzSpace, testSpace, entity, neighbour, tmpLocalMatrices[1], systemMatrix);
+ addToMatrix(ansatzSpace, testSpace, neighbour, entity, tmpLocalMatrices[2], systemMatrix);
+ addToMatrix(ansatzSpace, testSpace, neighbour, neighbour, tmpLocalMatrices[3], systemMatrix);
+ } // done visit only once
} // end if inner intersection
else if (!intersection.neighbor() && intersection.boundary()) // if boundary intersection
{
const unsigned int boundaryId = intersection.boundaryId();
- // if dirichlet boundary intersection
- if (boundaryId == 2) {
+ // // if dirichlet boundary intersection
+ // if( boundaryId == 2 )
+ // {
+
+ localDirichletOperator_.applyLocal(localAnsatzBaseFunctionSetEn,
+ localTesBaseFunctionSetEn,
+ intersection,
+ tmpLocalMatrices[0],
+ tmpLocalMatricesContainer[1]);
- } // end if dirichlet boundary intersection
- else if (boundaryId == 3) // if neumann boundary intersection
- {
+ addToMatrix(ansatzSpace, testSpace, entity, entity, tmpLocalMatrices[0], systemMatrix);
+
+ // } // end if dirichlet boundary intersection
+ // else if( boundaryId == 3 ) // if neumann boundary intersection
+ // {
+
+ // } // end if neumann boundary intersection
- } // end if neumann boundary intersection
} // end if boundary intersection
} // done loop over all intersections
- }
+ } // end method assembleLocal
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
+ template <class AnsatzSpaceType, class TestSpaceType, class EntityType, class LocalMatrixType, class SystemMatrixType>
+ void addToMatrix(const AnsatzSpaceType& ansatzSpace, const TestSpaceType& testSpace, const EntityType& ansatzEntity,
+ const EntityType& testEntity, const LocalMatrixType& localMatrix,
+ SystemMatrixType& systemMatrix) const
{
- unsigned int rows = ansatzSpace.baseFunctionSet().local(entity).size();
- unsigned int cols = testSpace.baseFunctionSet().local(entity).size();
+ unsigned int rows = ansatzSpace.baseFunctionSet().local(ansatzEntity).size();
+ unsigned int cols = testSpace.baseFunctionSet().local(testEntity).size();
for (unsigned int i = 0; i < rows; ++i) {
for (unsigned int j = 0; j < cols; ++j) {
- const unsigned int globalI = ansatzSpace.map().toGlobal(entity, i);
- const unsigned int globalJ = testSpace.map().toGlobal(neighbour, j);
+ const unsigned int globalI = ansatzSpace.map().toGlobal(ansatzEntity, i);
+ const unsigned int globalJ = testSpace.map().toGlobal(testEntity, j);
- matrix[globalI][globalJ] += localMatrix[i][j];
+ systemMatrix[globalI][globalJ] += localMatrix[i][j];
}
}
} // end method addToMatrix
- const InnerLocalOperatorType innerLocalOperator_;
- const DirichletLocalOperatorType dirichletLocalOperator_;
- const NeumannLocalOperatorType neumannLocalOperator_;
+ const LocalInnerOperatorType localInnerOperator_;
+ const LocalDirichletOperatorType localDirichletOperator_;
+ // const LocalNeumannOperatorType localNeumannOperator;
}; // end class Matrix
diff --git a/dune/detailed-discretizations/assembler/local/codim1/vector.hh b/dune/detailed-discretizations/assembler/local/codim1/vector.hh
new file mode 100644
index 0000000000000000000000000000000000000000..256975e1078cb01b5d92e7c74e8ff63789f3eb58
--- /dev/null
+++ b/dune/detailed-discretizations/assembler/local/codim1/vector.hh
@@ -0,0 +1,146 @@
+#ifndef DUNE_DETAILED_DISCRETIZATIONS_ASSEMLBER_LOCAL_CODIM1_VECTOR_HH
+#define DUNE_DETAILED_DISCRETIZATIONS_ASSEMLBER_LOCAL_CODIM1_VECTOR_HH
+
+// std includes
+#include <vector>
+
+namespace Dune {
+
+namespace DetailedDiscretizations {
+
+namespace Assembler {
+
+namespace Local {
+
+namespace Codim1 {
+
+/**
+ \todo Add neumann boundary treatment
+ \todo When adding neumann: think of numTmpObjectsRequired()!
+ \todo Add penalty parameter
+ **/
+template <class LocalFunctionalImp>
+class Vector
+{
+public:
+ typedef LocalFunctionalImp LocalFunctionalType;
+
+ typedef Vector<LocalFunctionalType> ThisType;
+
+ typedef typename LocalFunctionalType::RangeFieldType RangeFieldType;
+
+ //! constructor
+ Vector(const LocalFunctionalType localFunctional)
+ : localFunctional_(localFunctional)
+ {
+ }
+
+private:
+ //! copy constructor
+ Vector(const ThisType& other)
+ : localFunctional_(other.localFunctional())
+ {
+ }
+
+public:
+ const LocalFunctionalType& localFunctional() const
+ {
+ return localFunctional_;
+ }
+
+ /**
+ \todo Add neumann treatment here!
+ **/
+ std::vector<unsigned int> numTmpObjectsRequired() const
+ {
+ std::vector<unsigned int> ret(2, 0);
+ // we require 1 tmp vector in this local assembler
+ ret[0] = 1;
+ // the functional itself requires that much local matrices
+ ret[1] = localFunctional_.numTmpObjectsRequired();
+ return ret;
+ }
+
+ template <class TestSpaceType, class EntityType, class SystemVectorType, class LocalVectorType>
+ void assembleLocal(const TestSpaceType& testSpace, const EntityType& entity, SystemVectorType& systemVector,
+ std::vector<std::vector<LocalVectorType>>& tmpLocalVectorsContainer) const
+ {
+ // get the local basefunction set
+ typedef typename TestSpaceType::BaseFunctionSetType::LocalBaseFunctionSetType LocalTesBaseFunctionSetType;
+
+ const LocalTesBaseFunctionSetType localTestBaseFunctionSet = testSpace.baseFunctionSet().local(entity);
+
+ // check tmp local vectors
+ assert(tmpLocalVectorsContainer.size() > 1);
+ std::vector<LocalVectorType>& tmpLocalVectors = tmpLocalVectorsContainer[0];
+ if (tmpLocalVectors.size() < 1) {
+ tmpLocalVectors.resize(1, LocalVectorType(testSpace.map().maxLocalSize(), RangeFieldType(0.0)));
+ }
+
+ // some types
+ typedef typename TestSpaceType::GridPartType GridPartType;
+
+ typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
+
+ typedef typename IntersectionIteratorType::Intersection IntersectionType;
+
+ typedef typename IntersectionType::EntityPointer EntityPointerType;
+
+ const GridPartType& gridPart = testSpace.gridPart();
+
+ const IntersectionIteratorType lastIntersection = gridPart.iend(entity);
+
+ // do loop over all intersections
+ for (IntersectionIteratorType intIt = entity.ibegin(); intIt != lastIntersection; ++intIt) {
+ const IntersectionType& intersection = *intIt;
+
+ if (!intersection.neighbor() && intersection.boundary()) // if boundary intersection
+ {
+ const unsigned int boundaryId = intersection.boundaryId();
+
+ // // if dirichlet boundary intersection
+ // if( boundaryId == 2 )
+ // {
+ localFunctional_.applyLocal(
+ localTestBaseFunctionSet, intersection, tmpLocalVectors[0], tmpLocalVectorsContainer[1]);
+
+ // write local vector to global
+ addToVector(testSpace, entity, tmpLocalVectors[0], systemVector);
+
+ // } // end if dirichlet boundary intersection
+ // else if( boundaryId == 3 ) // if neumann boundary intersection
+ // {
+ // } // end if neumann boundary intersection
+ } // end if boundary intersection
+ } // done loop over all intersections
+ } // end method assembleLocal
+
+private:
+ //! assignment operator
+ ThisType& operator=(const ThisType&);
+
+ template <class TestSpaceType, class EntityType, class LocalVectorType, class SystemVectorType>
+ void addToVector(const TestSpaceType& testSpace, const EntityType& entity, const LocalVectorType& localVector,
+ SystemVectorType& systemVector) const
+ {
+ for (unsigned int j = 0; j < testSpace.baseFunctionSet().local(entity).size(); ++j) {
+ const unsigned int globalJ = testSpace.map().toGlobal(entity, j);
+
+ systemVector[globalJ] += localVector[j];
+ }
+ } // end method addToVector
+
+ const LocalFunctionalType localFunctional_;
+}; // end class Vector
+
+} // end namespace Codim1
+
+} // end namespace Local
+
+} // end namespace Assembler
+
+} // end namespace DetailedDiscretizations
+
+} // end namespace Dune
+
+#endif // DUNE_DETAILED_DISCRETIZATIONS_ASSEMLBER_LOCAL_CODIM1_VECTOR_HH
diff --git a/dune/detailed-discretizations/discretefunctional/local/Makefile.am b/dune/detailed-discretizations/discretefunctional/local/Makefile.am
index d931f05f8bee4396e467ddedc6e45a67a87fe4c1..827ea7f1590c698b75acb3956807b61eb58794e7 100644
--- a/dune/detailed-discretizations/discretefunctional/local/Makefile.am
+++ b/dune/detailed-discretizations/discretefunctional/local/Makefile.am
@@ -1,3 +1,3 @@
-SUBDIRS = codim0
+SUBDIRS = codim0 codim1
include $(top_srcdir)/am/global-rules
diff --git a/dune/detailed-discretizations/discretefunctional/local/codim1/Makefile.am b/dune/detailed-discretizations/discretefunctional/local/codim1/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..477d8bc45425b25abe527b7bb71e53decd5262b5
--- /dev/null
+++ b/dune/detailed-discretizations/discretefunctional/local/codim1/Makefile.am
@@ -0,0 +1,4 @@
+discretefunctional_local_codim1_includedir = $(includedir)/dune/detailed-discretizations/discretefunctional/local/codim1
+discretefunctional_local_codim1_include_HEADERS = integral.hh
+
+include $(top_srcdir)/am/global-rules
diff --git a/dune/detailed-discretizations/discretefunctional/local/codim1/integral.hh b/dune/detailed-discretizations/discretefunctional/local/codim1/integral.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3bcb3cbf319e35e229a714a1e28d828ccf9258e6
--- /dev/null
+++ b/dune/detailed-discretizations/discretefunctional/local/codim1/integral.hh
@@ -0,0 +1,259 @@
+/**
+ \file integral.hh
+ **/
+
+#ifndef DUNE_DETAILED_DISCRETIZATIONS_DISCRETEFUNCTIONAL_LOCAL_CODIM1_INTEGRAL_HH
+#define DUNE_DETAILED_DISCRETIZATIONS_DISCRETEFUNCTIONAL_LOCAL_CODIM1_INTEGRAL_HH
+
+// dune-common includes
+#include <dune/common/dynmatrix.hh>
+
+// dune fem includes
+#include <dune/fem/quadrature/cachingquadrature.hh>
+
+// dune-helper-tools includes
+#include <dune/helper-tools/common/vector.hh>
+
+namespace Dune {
+
+namespace DetailedDiscretizations {
+
+namespace DiscreteFunctional {
+
+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;
+
+ Integral(const LocalEvaluationType localEvaluation)
+ : localEvaluation_(localEvaluation)
+ {
+ }
+
+ //! copy constructor
+ Integral(const ThisType& other)
+ : localEvaluation_(other.localEvaluation())
+ {
+ }
+
+ LocalEvaluationType localEvaluation() const
+ {
+ return localEvaluation_;
+ }
+
+ unsigned int numTmpObjectsRequired() const
+ {
+ return 1;
+ }
+
+ template <class LocalTestBaseFunctionSetType, class IntersectionType, class LocalVectorType>
+ void applyLocal(const LocalTestBaseFunctionSetType localTestBaseFunctionSet, const IntersectionType& intersection,
+ LocalVectorType& localVector, std::vector<LocalVectorType>& tmpLocalVectors) const
+ {
+ // some types
+ typedef typename LocalTestBaseFunctionSetType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
+
+ typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
+
+ typedef Dune::CachingQuadrature<GridPartType, 0> VolumeQuadratureType;
+
+ // some stuff
+ const unsigned int size = localTestBaseFunctionSet.size();
+ const unsigned int quadratureOrder = localEvaluation_.order() + localTestBaseFunctionSet.order();
+ const VolumeQuadratureType volumeQuadrature(localTestBaseFunctionSet.entity(), quadratureOrder);
+ const unsigned int numberOfQuadraturePoints = volumeQuadrature.nop();
+
+ // make sure target vector is big enough
+ assert(localVector.size() >= size);
+
+ // clear target vector
+ Dune::HelperTools::Common::Vector::clear(localVector);
+
+ // check tmp local vectors
+ if (tmpLocalVectors.size() < 1) {
+ tmpLocalVectors.resize(1,
+ LocalVectorType(localTestBaseFunctionSet.baseFunctionSet().space().map().maxLocalSize(),
+ RangeFieldType(0.0)));
+ }
+
+ // do loop over all quadrature points
+ for (unsigned int q = 0; q < numberOfQuadraturePoints; ++q) {
+ // local coordinate
+ const DomainType x = volumeQuadrature.point(q);
+
+ // integration factors
+ const double integrationFactor = localTestBaseFunctionSet.entity().geometry().integrationElement(x);
+ const double quadratureWeight = volumeQuadrature.weight(q);
+
+ // evaluate the local evaluation
+ localEvaluation_.evaluateLocal(localTestBaseFunctionSet, x, tmpLocalVectors[0]);
+
+ // compute integral
+ for (unsigned int i = 0; i < size; ++i) {
+ localVector[i] += tmpLocalVectors[0][i] * integrationFactor * quadratureWeight;
+ }
+ } // done loop over all quadrature points
+ } // end method applyLocal
+
+private:
+ //! assignment operator
+ ThisType& operator=(const ThisType& other);
+
+ const LocalEvaluationType localEvaluation_;
+
+}; // end class Integral
+
+// template< class InducingOperatorImp, class InducingDiscreteFunctionImp >
+// class IntegralInduced
+//{
+// public:
+
+// typedef InducingOperatorImp
+// InducingOperatorType;
+
+// typedef InducingDiscreteFunctionImp
+// InducingDiscreteFunctionType;
+
+// typedef IntegralInduced< InducingOperatorType, InducingDiscreteFunctionType >
+// ThisType;
+
+// typedef typename InducingDiscreteFunctionType::RangeFieldType
+// RangeFieldType;
+
+// typedef typename InducingDiscreteFunctionType::DomainType
+// DomainType;
+
+// IntegralInduced( const InducingOperatorType& inducingOperator,
+// const InducingDiscreteFunctionType& inducingDiscreteFunction )
+// : inducingOperator_( inducingOperator ),
+// inducingDiscreteFunction_( inducingDiscreteFunction )
+// {
+// }
+
+// //! copy constructor
+// IntegralInduced( const ThisType& other )
+// : inducingOperator_( other.inducingOperator() ),
+// inducingDiscreteFunction_( other.inducingDiscreteFunction() )
+// {
+// }
+
+// InducingOperatorType inducingOperator() const
+// {
+// return inducingOperator_;
+// }
+
+// InducingDiscreteFunctionType inducingDiscreteFunction() const
+// {
+// return inducingDiscreteFunction_;
+// }
+
+// unsigned int numTmpObjectsRequired() const
+// {
+// return 0;
+// }
+
+// template< class LocalTestBaseFunctionSetType, class LocalVectorType >
+// void applyLocal( const LocalTestBaseFunctionSetType& localTestBaseFunctionSet,
+// LocalVectorType& localVector,
+// std::vector< LocalVectorType >& ) const
+// {
+// // some types
+// typedef typename LocalTestBaseFunctionSetType::DiscreteFunctionSpaceType
+// DiscreteFunctionSpaceType;
+
+// typedef typename DiscreteFunctionSpaceType::GridPartType
+// GridPartType;
+
+// typedef Dune::CachingQuadrature< GridPartType, 0 >
+// VolumeQuadratureType;
+
+// typedef typename LocalTestBaseFunctionSetType::EntityType
+// EntityType;
+
+// typedef typename InducingDiscreteFunctionType::ConstLocalFunctionType
+// InducingLocalFunctionType;
+
+// typedef typename Dune::DetailedDiscretizations::BaseFunctionSet::Local::Wrapper< InducingLocalFunctionType >
+// InducingBaseFunctionSetType;
+
+// typedef Dune::DynamicMatrix< RangeFieldType >
+// LocalMatrixType;
+
+// const EntityType& entity = localTestBaseFunctionSet.entity();
+
+// // wrap inducing local function
+// const InducingLocalFunctionType inducingLocalFunction = inducingDiscreteFunction_.localFunction( entity );
+// const InducingBaseFunctionSetType inducingBaseFunctionSet( inducingLocalFunction );
+
+// // some stuff
+// const unsigned int size = localTestBaseFunctionSet.size();
+// const unsigned int quadratureOrder =
+// inducingOperator_.localEvaluation().order() + inducingLocalFunction.order() + localTestBaseFunctionSet.order();
+// const VolumeQuadratureType volumeQuadrature( entity, quadratureOrder );
+// const unsigned int numberOfQuadraturePoints = volumeQuadrature.nop();
+
+// // make sure target vector is big enough
+// assert( localVector.size() >= size );
+
+// // clear target vector
+// Dune::HelperTools::Common::Vector::clear( localVector );
+
+// // do loop over all quadrature points
+// LocalMatrixType tmpMatrix( 1, size );
+// for( unsigned int q = 0; q < numberOfQuadraturePoints; ++q )
+// {
+// // local coordinate
+// const DomainType x = volumeQuadrature.point( q );
+
+// // integration factors
+// const double integrationFactor = entity.geometry().integrationElement( x );
+// const double quadratureWeight = volumeQuadrature.weight( q );
+
+// // evaluate the local evaluation
+// inducingOperator_.localEvaluation().evaluateLocal( inducingBaseFunctionSet,
+// localTestBaseFunctionSet,
+// x,
+// tmpMatrix );
+
+// // compute integral
+// for( unsigned int i = 0; i < size; ++i )
+// {
+// localVector[i] += tmpMatrix[0][i] * integrationFactor * quadratureWeight;
+// }
+// } // done loop over all quadrature points
+
+// } // end method applyLocal
+
+// private:
+// //! assignment operator
+// ThisType& operator=( const ThisType& );
+
+// const InducingOperatorType inducingOperator_;
+// const InducingDiscreteFunctionType inducingDiscreteFunction_;
+
+//}; // end class IntegralInduced
+
+} // end namespace Codim1
+
+} // end namespace Local
+
+} // end namespace DiscreteFunctional
+
+} // end namespace DetailedDiscretizations
+
+} // end namespace Dune
+
+#endif // end DUNE_DETAILED_DISCRETIZATIONS_DISCRETEFUNCTIONAL_LOCAL_CODIM0_INTEGRAL_HH
diff --git a/dune/detailed-discretizations/discreteoperator/local/Makefile.am b/dune/detailed-discretizations/discreteoperator/local/Makefile.am
index d931f05f8bee4396e467ddedc6e45a67a87fe4c1..827ea7f1590c698b75acb3956807b61eb58794e7 100644
--- a/dune/detailed-discretizations/discreteoperator/local/Makefile.am
+++ b/dune/detailed-discretizations/discreteoperator/local/Makefile.am
@@ -1,3 +1,3 @@
-SUBDIRS = codim0
+SUBDIRS = codim0 codim1
include $(top_srcdir)/am/global-rules
diff --git a/dune/detailed-discretizations/discreteoperator/local/codim1/Makefile.am b/dune/detailed-discretizations/discreteoperator/local/codim1/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..bf452c3455052a0a9b43522daf4bc33f046e5cfa
--- /dev/null
+++ b/dune/detailed-discretizations/discreteoperator/local/codim1/Makefile.am
@@ -0,0 +1,4 @@
+discreteoperator_local_codim1_includedir = $(includedir)/dune/detailed-discretizations/discreteoperator/local/codim1
+discreteoperator_local_codim1_include_HEADERS = boundaryintegral.hh innerintegral.hh
+
+include $(top_srcdir)/am/global-rules
diff --git a/dune/detailed-discretizations/discreteoperator/local/codim1/boundaryintegral.hh b/dune/detailed-discretizations/discreteoperator/local/codim1/boundaryintegral.hh
index b6a99ace767f9c2eefa4090fb38e2964437d7c85..7ed0c38bbf0b217ab75d9e8551922ff046ec32d4 100644
--- a/dune/detailed-discretizations/discreteoperator/local/codim1/boundaryintegral.hh
+++ b/dune/detailed-discretizations/discreteoperator/local/codim1/boundaryintegral.hh
@@ -80,9 +80,6 @@ public:
return 1;
}
- /**
- \todo Rename Entity -> En, Neighbour -> Ne
- **/
template <class LocalAnsatzBaseFunctionSetType, class LocalTestBaseFunctionSetType, class IntersectionType,
class LocalMatrixType>
void applyLocal(const LocalAnsatzBaseFunctionSetType& localAnsatzBaseFunctionSet,
@@ -131,13 +128,12 @@ public:
// evaluate the local operation
localEvaluation_.evaluate(
- localAnsatzBaseFunctionSet, localTestBaseFunctionSet, intersection, x, tmpLocalMatrices[0]); /*NeNe*/
+ localAnsatzBaseFunctionSet, localTestBaseFunctionSet, intersection, x, tmpLocalMatrices[0]);
// compute integral (see below)
addToIntegral(tmpLocalMatrices[0], integrationFactor, quadratureWeight, rows, cols, localMatrix);
} // done loop over all quadrature points
-
} // end method applyLocal
private:
diff --git a/dune/detailed-discretizations/evaluation/local/binary/Makefile.am b/dune/detailed-discretizations/evaluation/local/binary/Makefile.am
index 4ba625dae47f36901edfcc238dd3890991a78df9..acdbff6c03f308e925538344ca9c7a071449def8 100644
--- a/dune/detailed-discretizations/evaluation/local/binary/Makefile.am
+++ b/dune/detailed-discretizations/evaluation/local/binary/Makefile.am
@@ -1,4 +1,4 @@
-evaluation_binary_includedir = $(includedir)/dune/detailed-discretizations/evaluation/binary
-evaluation_binary_include_HEADERS = elliptic.hh ipdgfluxes.hh
+evaluation_local_binary_includedir = $(includedir)/dune/detailed-discretizations/evaluation/local/binary
+evaluation_local_binary_include_HEADERS = elliptic.hh ipdgfluxes.hh
include $(top_srcdir)/am/global-rules
diff --git a/dune/detailed-discretizations/evaluation/local/quaternary/Makefile.am b/dune/detailed-discretizations/evaluation/local/quaternary/Makefile.am
index 70f5521b9676b4d744b975f28c00c725c30bf5a2..46aac1b97ac966ce31a98d29361ef1fb3940050f 100644
--- a/dune/detailed-discretizations/evaluation/local/quaternary/Makefile.am
+++ b/dune/detailed-discretizations/evaluation/local/quaternary/Makefile.am
@@ -1,4 +1,4 @@
-evaluation_quaternary_includedir = $(includedir)/dune/detailed-discretizations/evaluation/quaternary
-evaluation_quaternary_include_HEADERS = ipdgfluxes.hh
+evaluation_local_quaternary_includedir = $(includedir)/dune/detailed-discretizations/evaluation/local/quaternary
+evaluation_local_quaternary_include_HEADERS = ipdgfluxes.hh
include $(top_srcdir)/am/global-rules
diff --git a/dune/detailed-discretizations/evaluation/local/unary/Makefile.am b/dune/detailed-discretizations/evaluation/local/unary/Makefile.am
index 32fe97af4fb9f6d66af800611d805b51e8c517d7..eee81ca982b0a6f9004dd7e6dad864f097c5346b 100644
--- a/dune/detailed-discretizations/evaluation/local/unary/Makefile.am
+++ b/dune/detailed-discretizations/evaluation/local/unary/Makefile.am
@@ -1,4 +1,4 @@
-evaluation_unary_includedir = $(includedir)/dune/detailed-discretizations/evaluation/unary
-evaluation_unary_include_HEADERS = scale.hh
+evaluation_local_unary_includedir = $(includedir)/dune/detailed-discretizations/evaluation/local/unary
+evaluation_local_unary_include_HEADERS = scale.hh
include $(top_srcdir)/am/global-rules
diff --git a/dune/detailed-discretizations/evaluation/local/unary/ipdgfluxes.hh b/dune/detailed-discretizations/evaluation/local/unary/ipdgfluxes.hh
new file mode 100644
index 0000000000000000000000000000000000000000..744c23c49ad92f684ec6c7e85149736571a9769e
--- /dev/null
+++ b/dune/detailed-discretizations/evaluation/local/unary/ipdgfluxes.hh
@@ -0,0 +1,137 @@
+#ifndef DUNE_DETAILED_DISCRETIZATIONS_EVALUATION_LOCAL_UNARY_IPDGFLUXES_HH
+#define DUNE_DETAILED_DISCRETIZATIONS_EVALUATION_LOCAL_UNARY_IPDGFLUXES_HH
+
+// dune-helper-tools includes
+#include <dune/helper-tools/function/runtime.hh>
+
+namespace Dune {
+
+namespace DetailedDiscretizations {
+
+namespace Evaluation {
+
+namespace Local {
+
+namespace Unary {
+
+/**
+ \todo Should be parameterized with the inducing and the dirichlet boundary function.
+ \todo Add neumann.
+ \todo Add penalty parameter.
+ \tparam FunctionSpaceImp
+ Type of the function space, where \f$f\f$ and \f$v\f$ live in.
+ **/
+template <class FunctionSpaceImp>
+class IPDGFluxes
+{
+public:
+ typedef FunctionSpaceImp FunctionSpaceType;
+
+ typedef IPDGFluxes<FunctionSpaceType> ThisType;
+
+ typedef typename FunctionSpaceType::DomainType DomainType;
+
+ typedef typename FunctionSpaceType::RangeType RangeType;
+
+ typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
+
+ typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
+
+ typedef Dune::HelperTools::Function::Runtime<FunctionSpaceType> InducingFunctionType;
+
+ typedef Dune::HelperTools::Function::Runtime<FunctionSpaceType> DirichletFunctionType;
+
+ //! constructor, takes the inducing and the dirichlet functions expressions as runtime parameters
+ IPDGFluxes(const std::string inducingExpression = "[1.0;1.0;1.0]",
+ const std::string dirichletExpression = "[0.0;0.0;0.0]", const int order = 1)
+ : inducingFunction_(inducingExpression)
+ , dirichletFunction_(dirichletExpression)
+ , order_(std::max(0, order))
+ {
+ }
+
+ //! copy constructor
+ IPDGFluxes(const ThisType& other)
+ : inducingFunction_(other.inducingFunction())
+ , dirichletFunction_(other.dirichletFunction())
+ , order_(other.order())
+ {
+ }
+
+ //! returns the inducing function
+ InducingFunctionType inducingFunction() const
+ {
+ return inducingFunction_;
+ }
+
+ DirichletFunctionType dirichletFunction() const
+ {
+ return dirichletFunction_;
+ }
+
+ unsigned int order() const
+ {
+ return order_;
+ }
+
+ /**
+ \attention Assumes ret to be cleared, since we do multiple +=
+ **/
+ template <class LocalTestBaseFunctionSetType, class IntersectionType, class LocalVectorType>
+ void evaluateLocal(const LocalTestBaseFunctionSetType& localTestBaseFunctionSet, const IntersectionType& intersection,
+ const DomainType& localPoint, LocalVectorType& ret) const
+ {
+ // some stuff
+ const DomainType globalPoint = intersection.geometry().global(localPoint);
+ const DomainType unitOuterNormal = intersection.unitOuterNormal();
+
+ // evaluate inducing and dirichlet function
+ RangeType inducingFunctionEvaluation(0.0);
+ inducingFunction_.evaluate(globalPoint, inducingFunctionEvaluation);
+ RangeType dirichletFunctionEvaluation(0.0);
+ inducingFunction_.evaluate(globalPoint, dirichletFunctionEvaluation);
+
+ // evaluate test basefunctionset
+ const unsigned int size = localTestBaseFunctionSet.size();
+ std::vector<RangeType> localTestBaseFunctionSetEvaluations(size, RangeType(0.0));
+ localTestBaseFunctionSet.evaluate(localPoint, localTestBaseFunctionSetEvaluations);
+ std::vector<JacobianRangeType> localTestBaseFunctionSetGradients(size, JacobianRangeType(0.0));
+ localTestBaseFunctionSet.jacobian(localPoint, localTestBaseFunctionSetGradients);
+
+ // evaluate penalty parameter
+ const RangeFieldType penaltyParameter = 1.0 / std::pow(intersection.geometry().volume(), 1.0);
+
+ // do loop over all basis functions
+ assert(ret.size() >= size);
+ for (unsigned int i = 0; i < size; ++i) {
+ {
+ const RangeType evaluationTimesDirichlet = localTestBaseFunctionSetEvaluations[i] * dirichletFunctionEvaluation;
+ ret[i] = penaltyParameter * evaluationTimesDirichlet;
+ }
+ {
+ const RangeFieldType gradientTimesNormal = localTestBaseFunctionSetGradients[i][0] * unitOuterNormal;
+ ret[i] += -1.0 * inducingFunctionEvaluation * gradientTimesNormal * dirichletFunctionEvaluation;
+ }
+ }
+ } // end method evaluateLocal
+
+private:
+ //! assignment operator
+ ThisType& operator=(const ThisType&);
+
+ const InducingFunctionType inducingFunction_;
+ const DirichletFunctionType dirichletFunction_;
+ const unsigned int order_;
+}; // end class IPDGFluxes
+
+} // end namespace Unary
+
+} // end namespace Local
+
+} // end namespace Evaluation
+
+} // end namespace DetailedDiscretizations
+
+} // end namespace Dune
+
+#endif // DUNE_DETAILED_DISCRETIZATIONS_EVALUATION_LOCAL_UNARY_IPDGFLUXES_HH
diff --git a/examples/elliptic_discontinuous_galerkin/main.cc b/examples/elliptic_discontinuous_galerkin/main.cc
index 82f61d14eae45b08cd39cd7025393c18e7a92dcf..4eeb21ace7f0d682ed181cc9867476c5ca1dc9e4 100644
--- a/examples/elliptic_discontinuous_galerkin/main.cc
+++ b/examples/elliptic_discontinuous_galerkin/main.cc
@@ -37,10 +37,13 @@
// dune-detailed-discretizations includes
#include <dune/detailed-discretizations/discretefunctionspace/continuous/lagrange.hh>
#include <dune/detailed-discretizations/evaluation/local/quaternary/ipdgfluxes.hh>
+#include <dune/detailed-discretizations/evaluation/local/unary/ipdgfluxes.hh>
#include <dune/detailed-discretizations/discreteoperator/local/codim1/innerintegral.hh>
#include <dune/detailed-discretizations/discreteoperator/local/codim1/boundaryintegral.hh>
+#include <dune/detailed-discretizations/discretefunctional/local/codim1/integral.hh>
#include <dune/detailed-discretizations/container/factory.hh>
-//#include <dune/detailed-discretizations/assembler/local/codim1/matrix.hh>
+#include <dune/detailed-discretizations/assembler/local/codim1/matrix.hh>
+#include <dune/detailed-discretizations/assembler/local/codim1/vector.hh>
// only ever do this in a .cc!
using namespace Dune::DetailedDiscretizations;
@@ -109,6 +112,10 @@ int main(int argc, char** argv)
const DirichletFacesIPDGEvaluationType dirichletFacesIPDGEvaluation("[1.0;1.0;1.0]");
+ typedef Dune::DetailedDiscretizations::Evaluation::Local::Unary::IPDGFluxes<FunctionSpaceType> IPDGEvaluationType;
+
+ const IPDGEvaluationType ipdgEvaluation("[1.0;1.0;1.0]", "[0.0;0.0;0.0]", 1);
+
// operator and functional
// typedef DiscreteOperator::Local::Codim0::Integral< EllipticEvaluationType >
@@ -116,11 +123,6 @@ int main(int argc, char** argv)
// const LocalEllipticOperatorType localEllipticOperator( ellipticEvaluation );
- // typedef DiscreteFunctional::Local::Codim0::Integral< ProductEvaluationType >
- // LocalL2FunctionalType;
-
- // const LocalL2FunctionalType localL2Functional( productEvaluation );
-
typedef DiscreteOperator::Local::Codim1::InnerIntegral<InnerFacesIPDGEvaluationType>
LocalIPDGInnerFacesOperatorType;
@@ -131,6 +133,15 @@ int main(int argc, char** argv)
const LocalIPDGDirichletFacesOperatorType localIPDGDirichletFacesOperator(dirichletFacesIPDGEvaluation);
+ // typedef DiscreteFunctional::Local::Codim0::Integral< ProductEvaluationType >
+ // LocalL2FunctionalType;
+
+ // const LocalL2FunctionalType localL2Functional( productEvaluation );
+
+ typedef DiscreteFunctional::Local::Codim1::Integral<IPDGEvaluationType> LocalIPDGFunctionalType;
+
+ const LocalIPDGFunctionalType localIPDGFunctional(ipdgEvaluation);
+
// matrix, rhs and solution storage
typedef Container::Matrix::Defaults<RangeFieldType, dimRange, dimRange>::BCRSMatrix MatrixFactory;
@@ -152,14 +163,24 @@ int main(int argc, char** argv)
// assembler
// typedef Assembler::Local::Codim0::Matrix< LocalEllipticOperatorType >
- // LocalMatrixAssemblerType;
+ // LocalCodim0MatrixAssemblerType;
- // const LocalMatrixAssemblerType localMatrixAssembler( localEllipticOperator );
+ // const LocalCodim0MatrixAssemblerType localCodim0MatrixAssembler( localEllipticOperator );
+
+ typedef Assembler::Local::Codim1::Matrix<LocalIPDGInnerFacesOperatorType, LocalIPDGDirichletFacesOperatorType>
+ LocalCodim1MatrixAssemblerType;
+
+ const LocalCodim1MatrixAssemblerType localCodim1MatrixAssembler(localIPDGInnerFacesOperator,
+ localIPDGDirichletFacesOperator);
// typedef Assembler::Local::Codim0::Vector< LocalL2FunctionalType >
- // LocalVectorAssemblerType;
+ // LocalCodim0VectorAssemblerType;
+
+ // const LocalCodim0VectorAssemblerType localCodim0VectorAssembler( localL2Functional );
+
+ typedef Assembler::Local::Codim1::Vector<LocalIPDGFunctionalType> LocalCodim1VectorAssemblerType;
- // const LocalVectorAssemblerType localVectorAssembler( localL2Functional );
+ const LocalCodim1VectorAssemblerType localCodim1VectorAssembler(localIPDGFunctional);
// typedef Assembler::System::Unconstrained< DiscreteH1Type, DiscreteH1Type >
// SystemAssemblerType;