Skip to content
Snippets Groups Projects
Commit 66a8b073 authored by Patrick Jaap's avatar Patrick Jaap
Browse files

Assembler for global basis transfer matrices

The method assembles the transfer matrx between fine and coarse grids.
Assumptions are
 - the grids are related but not necessarily direct descendants
 - the global dune-functions basis are scalar, i.e., leaf nodes in the basis tree

For evaluating local functions on other geometries a small wrapper is written.
For chaining sequencial father geometries another small wrapper is written.

Two unit tests:
- first, test wether the transfermatrix interpolates coefficients from corase to fine
  for lagrange basis
- second, test wether a jump over two grid levels can be reproduced by two simple transfers
parent f0f5e72e
No related branches found
No related tags found
1 merge request!40Assembler for global basis transfer matrices
Pipeline #12310 passed
......@@ -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
......@@ -3,6 +3,7 @@ set(GRID_BASED_TESTS
basisgridfunctiontest
basisinterpolatortest
basisinterpolationmatrixtest
assembletransferoperatortest
boundarypatchtest
boundarypatchprolongatortest
coarsegridfunctionwrappertest
......
#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;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment