diff --git a/dune/fufem/assemblers/basisinterpolationmatrixassembler.hh b/dune/fufem/assemblers/basisinterpolationmatrixassembler.hh
index 91c2bd9b8fd424caef4d85c2ceb26c31a940d75e..8478c9b6a9a5a91c2c9723677b7eac43686b954b 100644
--- a/dune/fufem/assemblers/basisinterpolationmatrixassembler.hh
+++ b/dune/fufem/assemblers/basisinterpolationmatrixassembler.hh
@@ -4,6 +4,8 @@
#define DUNE_FUFEM_ASSEMBLERS_BASIS_INTERPOLATION_MATRIX_ASSEMBLER_HH
#include <type_traits>
+#include <typeinfo>
+
#include <dune/common/bitsetvector.hh>
#include <dune/istl/matrixindexset.hh>
@@ -215,4 +217,281 @@ static void assembleBasisInterpolationMatrix(MatrixType& matrix, const BasisType
}
}
+
+
+
+
+
+/** \brief Wrapper for evaluating local functions in another geometry */
+template<class Function, class Geometry>
+class OtherGeometryWrapper
+{
+
+public:
+
+ using Traits = typename Function::Traits;
+
+ OtherGeometryWrapper(const Function& function, const Geometry& geometry)
+ : function_(function)
+ , geometry_(geometry)
+ {}
+
+ // forward the evaluation into the other geometry
+ template<class X, class Y>
+ void evaluate(const X& x, Y& y) const
+ {
+ function_.evaluate(geometry_.global(x), y);
+ }
+
+private:
+
+ const Function function_;
+ const Geometry geometry_;
+};
+
+
+/** \brief Wrapper chaining sequencial father geometries */
+template<class GeometryContainer>
+class GeometryInFathers
+{
+
+public:
+
+ GeometryInFathers(const GeometryContainer& geometryContainer)
+ : geometryContainer_(geometryContainer)
+ {}
+
+ // forward the global() method to all single geometries
+ template<class X>
+ auto global(const X& x) const
+ {
+ auto y = x;
+ for (const auto& g : geometryContainer_)
+ {
+ y = g.global(y);
+ }
+ return y;
+ }
+
+private:
+
+ const GeometryContainer geometryContainer_;
+};
+
+
+
+/** \brief Assemble the transfer matrix for multigrid methods
+ *
+ * The coarse and fine basis has to be related: in a way that the fine grid is descendant of the
+ * coarse grid (not necessarily direct descandant).
+ * The matrix is assumed to be in block structure with correct block size.
+ * The two basis are assumed to be scalar dune-functions basis (i.e. leaf nodes in the basis tree).
+ *
+ * \param [out] matrix the block matrix corresponding to the basis transfer operator
+ * \param [in] coarseBasis scalar global basis on the coarse grid
+ * \param [in] fineasis scalar global basis on the fine grid
+ * \param [in] tolerance (optional) interpolation threshold. Default is 1e-12.
+ */
+template<class MatrixType, class CoarseBasis, class FineBasis>
+static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
+ const CoarseBasis& coarseBasis,
+ const FineBasis& fineBasis,
+ const double tolerance = 1e-12)
+{
+ const auto& coarseGridView = coarseBasis.gridView();
+ const auto& fineGridView = fineBasis.gridView();
+
+ auto coarseLocalView = coarseBasis.localView();
+ auto fineLocalView = fineBasis.localView();
+
+ const auto rows = fineBasis.size();
+ const auto cols = coarseBasis.size();
+
+ if ( not fineLocalView.tree().isLeaf or not coarseLocalView.tree().isLeaf )
+ DUNE_THROW(Dune::Exception, "currently, we can handle power nodes only");
+
+ using CoarseFEType = typename decltype(coarseLocalView)::Tree::FiniteElement;
+ using FunctionBaseClass = typename Dune::LocalFiniteElementFunctionBase<CoarseFEType>::type;
+ using LocalBasisWrapper = LocalBasisComponentWrapper<typename CoarseFEType::Traits::LocalBasisType, FunctionBaseClass>;
+
+ matrix.setSize(rows,cols);
+ matrix.setBuildMode(MatrixType::random);
+
+
+ // ///////////////////////////////////////////
+ // Determine the indices present in the matrix
+ // ///////////////////////////////////////////
+
+ // Only handle every dof once
+ Dune::BitSetVector<1> processed(fineBasis.size());
+
+ Dune::MatrixIndexSet indices(rows, cols);
+
+ // loop over all fine grid elements
+ for ( const auto& fineElement : elements(fineGridView) )
+ {
+ // find a coarse element that is the parent of the fine element
+ // start in the fine grid and track the geometry mapping
+ auto fE = fineElement;
+ // collect all geometryInFather's on the way
+ std::vector<decltype(fineElement.geometryInFather())> geometryInFathersVector;
+ while ( not coarseGridView.contains(fE) and fE.hasFather() )
+ {
+ // add the geometry to the container
+ geometryInFathersVector.push_back(fE.geometryInFather());
+ // step up one level
+ fE = fE.father();
+ }
+
+ // did we really end up in coarseElement?
+ if ( not coarseGridView.contains(fE) )
+ DUNE_THROW( Dune::Exception, "There is a fine element without a parent in the coarse grid!");
+
+ const auto& coarseElement = fE;
+
+ const auto geometryInFathers = GeometryInFathers<decltype(geometryInFathersVector)>(geometryInFathersVector);
+
+ coarseLocalView.bind(coarseElement);
+ fineLocalView.bind(fineElement);
+
+ // get the root as reference
+ const auto& node = fineLocalView.tree();
+
+ // get the local basis
+ const auto& localBasis = node.finiteElement().localBasis();
+ using Range = typename std::decay_t<decltype(localBasis)>::Traits::RangeType;
+
+ std::vector<Range> values(localBasis.size());
+
+ LocalBasisWrapper basisFctj(node.finiteElement().localBasis(),0);
+
+ const auto& coarseLocalBasis = coarseLocalView.tree().finiteElement().localBasis();
+
+ for (size_t j=0; j<coarseLocalBasis.size(); j++)
+ {
+ // wrap each local basis function as a local function.
+ basisFctj.setIndex(j);
+
+ // transform the local fine function to a local coarse function
+ const auto coarseBaseFctj = OtherGeometryWrapper<decltype(basisFctj),decltype(geometryInFathers)>(basisFctj, geometryInFathers);
+
+ // Interpolate j^th base function by the fine basis
+ node.finiteElement().localInterpolation().interpolate(coarseBaseFctj, values);
+
+ // get the matrix col index
+ const auto& globalCoarseIndex = coarseLocalView.index(j);
+
+ // set the matrix indices
+ for (size_t i=0; i<localBasis.size(); i++) {
+
+ if ( std::abs(values[i]) < tolerance )
+ continue;
+
+ // get the matrix row index
+ const auto& globalFineIndex = fineLocalView.index(i);
+
+ if (processed[globalFineIndex][0])
+ continue;
+
+ indices.add(globalFineIndex, globalCoarseIndex);
+ }
+ }
+
+ // now the all dof's of this fine element to 'processed'
+ for (size_t i=0; i<localBasis.size(); i++)
+ {
+ const auto& globalFineIndex = fineLocalView.index(i);
+ processed[globalFineIndex] = true;
+ }
+ }
+
+ indices.exportIdx(matrix);
+ matrix = 0;
+
+ // reset all dof's
+ processed.unsetAll();
+
+ //////////////////////
+ // Now set the values
+ //////////////////////
+
+ for ( const auto& fineElement : elements(fineGridView) )
+ {
+ // find a coarse element that is the parent of the fine element
+ // start in the fine grid and track the geometry mapping
+ auto fE = fineElement;
+ // collect all geometryInFather's on the way
+ std::vector<decltype(fineElement.geometryInFather())> geometryInFathersVector;
+ while ( not coarseGridView.contains(fE) and fE.hasFather() )
+ {
+ // add the geometry to the container
+ geometryInFathersVector.push_back(fE.geometryInFather());
+ // step up one level
+ fE = fE.father();
+ }
+
+ // did we really end up in coarseElement?
+ if ( not coarseGridView.contains(fE) )
+ DUNE_THROW( Dune::Exception, "There is a fine element without a parent in the coarse grid!");
+
+ const auto& coarseElement = fE;
+
+ const auto geometryInFathers = GeometryInFathers<decltype(geometryInFathersVector)>(geometryInFathersVector);
+
+ coarseLocalView.bind(coarseElement);
+ fineLocalView.bind(fineElement);
+
+ // get the root as reference
+ const auto& node = fineLocalView.tree();
+
+ // get the local basis
+ const auto& localBasis = node.finiteElement().localBasis();
+ using Range = typename std::decay_t<decltype(localBasis)>::Traits::RangeType;
+
+ std::vector<Range> values(localBasis.size());
+
+ LocalBasisWrapper basisFctj(node.finiteElement().localBasis(),0);
+
+ const auto& coarseLocalBasis = coarseLocalView.tree().finiteElement().localBasis();
+
+ for (size_t j=0; j<coarseLocalBasis.size(); j++)
+ {
+ // wrap each local basis function as a local function.
+ basisFctj.setIndex(j);
+
+ // transform the local fine function to a local coarse function
+ const auto coarseBaseFctj = OtherGeometryWrapper<decltype(basisFctj),decltype(geometryInFathers)>(basisFctj, geometryInFathers);
+
+ // Interpolate j^th base function by the fine basis
+ node.finiteElement().localInterpolation().interpolate(coarseBaseFctj, values);
+
+ // get the matrix col index
+ const auto& globalCoarseIndex = coarseLocalView.index(j);
+
+ // set the matrix indices
+ for (size_t i=0; i<localBasis.size(); i++) {
+
+ if ( std::abs(values[i]) < tolerance )
+ continue;
+
+ // get the matrix row index
+ const auto& globalFineIndex = fineLocalView.index(i);
+
+ if (processed[globalFineIndex][0])
+ continue;
+
+ Dune::MatrixVector::addToDiagonal(matrix[globalFineIndex][globalCoarseIndex],values[i]);
+ }
+ }
+
+ // now the all dof's of this fine element to 'processed'
+ for (size_t i=0; i<localBasis.size(); i++)
+ {
+ const auto& globalFineIndex = fineLocalView.index(i);
+ processed[globalFineIndex] = true;
+ }
+ }
+}
+
+
#endif
diff --git a/dune/fufem/test/CMakeLists.txt b/dune/fufem/test/CMakeLists.txt
index a45572682ad4aa17372c51d9e1e59c7488aae60e..2c6685b61159a752c2e51a278403fd4e4ad37ba1 100644
--- a/dune/fufem/test/CMakeLists.txt
+++ b/dune/fufem/test/CMakeLists.txt
@@ -3,6 +3,7 @@ set(GRID_BASED_TESTS
basisgridfunctiontest
basisinterpolatortest
basisinterpolationmatrixtest
+ assembletransferoperatortest
boundarypatchtest
boundarypatchprolongatortest
coarsegridfunctionwrappertest
diff --git a/dune/fufem/test/assembletransferoperatortest.cc b/dune/fufem/test/assembletransferoperatortest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..bbad6a37e9b4e193338ec17f8e664372ec74c5db
--- /dev/null
+++ b/dune/fufem/test/assembletransferoperatortest.cc
@@ -0,0 +1,162 @@
+#include <config.h>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/test/testsuite.hh>
+#include <dune/istl/bcrsmatrix.hh>
+
+#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
+#include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh>
+
+#include "common.hh"
+
+using namespace Dune;
+
+template<int order, class GW0, class GW1>
+bool testLagrange(const GW0& gw0, const GW1& gw1)
+{
+ // our test funcion for interpolation
+ auto f = [](const auto& x)
+ {
+ // const function for order 0
+ if ( order == 0 )
+ return 1.;
+
+ // for the other cases affine function
+ auto y = 1.;
+ for( auto xs : x )
+ y += xs;
+ return y;
+ };
+
+ auto coarseLagrange = Functions::BasisFactory::makeBasis(gw0, Functions::BasisFactory::lagrange<order>());
+ auto fineLagrange = Functions::BasisFactory::makeBasis(gw1, Functions::BasisFactory::lagrange<order>());
+
+ BCRSMatrix<FieldMatrix<double,1,1>> matrix;
+ BlockVector<FieldVector<double,1>> coarseCoeffs, fineCoeffs;
+
+ // assemble the matrix
+ assembleGlobalBasisTransferMatrix(matrix,coarseLagrange,fineLagrange);
+
+ // interpolate f on the coarse grid
+ Functions::interpolate(coarseLagrange,coarseCoeffs,f);
+
+ // interpolate f on the fine grid
+ Functions::interpolate(fineLagrange,fineCoeffs,f);
+
+ // now we should have
+ // fineCoeffs - matrix * coarseCoeffs == 0
+ matrix.mmv(coarseCoeffs,fineCoeffs);
+
+ std::cout << " test lagrage " << GW0::dimension << "D with order " << order << " : error = " << fineCoeffs.two_norm2()<< std::endl;
+
+ return fineCoeffs.two_norm2() < 1e-12;
+}
+
+
+template<int order, class GW0, class GW1, class GW2>
+bool testTwoLevelLagrange(const GW0& gw0, const GW1& gw1, const GW2& gw2)
+{
+
+ // create 3 global basis
+ auto coarseLagrange = Functions::BasisFactory::makeBasis(gw0, Functions::BasisFactory::lagrange<order>());
+ auto semiLagrange = Functions::BasisFactory::makeBasis(gw1, Functions::BasisFactory::lagrange<order>());
+ auto fineLagrange = Functions::BasisFactory::makeBasis(gw2, Functions::BasisFactory::lagrange<order>());
+
+ // matrices for the conbined jump and two simple jumps
+ BCRSMatrix<FieldMatrix<double,1,1>> matrixTwoJumps, matrixJumpOne, matrixJumpTwo;
+
+ // assemble the matrices
+ assembleGlobalBasisTransferMatrix(matrixTwoJumps,coarseLagrange,fineLagrange);
+ assembleGlobalBasisTransferMatrix(matrixJumpOne,coarseLagrange,semiLagrange);
+ assembleGlobalBasisTransferMatrix(matrixJumpTwo,semiLagrange,fineLagrange);
+
+ // now check wether the product of two jumps is the conbined jump
+ // compute matrixTwoJumps -= matrixJumpTwo*matrixJumpOne
+ for (auto row0 = matrixJumpTwo.begin(); row0 != matrixJumpTwo.end(); row0++)
+ for (auto entry0 = row0->begin(); entry0 != row0->end(); entry0++)
+ for (auto row1 = matrixJumpOne.begin(); row1 != matrixJumpOne.end(); row1++)
+ for(auto entry1 = row1->begin(); entry1 != row1->end(); entry1++)
+ if (row1.index() == entry0.index() )
+ matrixTwoJumps[row0.index()][entry1.index()] -= (*entry0)*(*entry1);
+
+
+ std::cout << " test 2-jump " << GW0::dimension << "D with order " << order << " : error = " << matrixTwoJumps.frobenius_norm2()<< std::endl;
+
+ return matrixTwoJumps.frobenius_norm2() < 1e-15;
+
+}
+
+
+
+
+// this tests the interpolation porperty of the matrix
+template<class GW0, class GW1>
+bool checkTransferMatrix(const GW0& gw0, const GW1& gw1)
+{
+ bool passed = true;
+
+ // test oder 0-3
+ passed = passed and testLagrange<0>(gw0,gw1);
+ passed = passed and testLagrange<1>(gw0,gw1);
+ passed = passed and testLagrange<2>(gw0,gw1);
+ // order 3 only for 2D grids
+ if ( GW0::dimension < 3 )
+ passed = passed and testLagrange<3>(gw0,gw1);
+
+ return passed;
+}
+
+
+// this test checks a jump over two levels
+template<class GW0, class GW1, class GW2>
+bool checkTwoLevelJump(const GW0& gw0, const GW1& gw1, const GW2& gw2)
+{
+ bool passed = true;
+
+ // test oder 0-3
+ passed = passed and testTwoLevelLagrange<0>(gw0,gw1,gw2);
+ passed = passed and testTwoLevelLagrange<1>(gw0,gw1,gw2);
+ passed = passed and testTwoLevelLagrange<2>(gw0,gw1,gw2);
+ // order 3 only for 2D grids
+ if ( GW0::dimension < 3 )
+ passed = passed and testTwoLevelLagrange<3>(gw0,gw1,gw2);
+
+ return passed;
+}
+
+
+struct Suite
+{
+ template<class GridType>
+ bool check(const GridType& grid)
+ {
+ bool passed = true;
+
+ auto maxLevel = grid.maxLevel();
+
+ auto fineGridView = grid.levelGridView(maxLevel);
+ auto semiGridView = grid.levelGridView(maxLevel-1);
+ auto coarseGridView = grid.levelGridView(maxLevel-2);
+
+ passed = passed and checkTransferMatrix(semiGridView,fineGridView);
+ passed = passed and checkTwoLevelJump(coarseGridView,semiGridView,fineGridView);
+
+ return passed;
+ }
+};
+
+int main(int argc, char** argv)
+{
+ Dune::MPIHelper::instance(argc, argv);
+
+ Suite tests;
+
+ bool passed = checkWithStructuredGrid<2>(tests, 3);
+ passed = passed and checkWithStructuredGrid<3>(tests, 3);
+
+ return passed ? 0 : 1;
+}