Merge branch 'feature/transfer-matrix' into 'master'

Assembler for global basis transfer matrices

See merge request !40

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

dune/fufem/test/CMakeLists.txt

@@ -3,6 +3,7 @@ set(GRID_BASED_TESTS
   basisgridfunctiontest
   basisinterpolatortest
   basisinterpolationmatrixtest
+  assembletransferoperatortest
   boundarypatchtest
   boundarypatchprolongatortest
   coarsegridfunctionwrappertest

dune/fufem/test/assembletransferoperatortest.cc

+#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;
+}