Commit 6d90eabe authored by Patrick Jaap's avatar Patrick Jaap

Merge branch 'feature/assemble-transfer-operator-with-powerbasis' into 'master'

assembleGlobalBasisTransferMatrix: Extend to dune-functions bases

See merge request !65
parents fc853df5 3d52dbb1
Pipeline #28193 passed with stage
in 13 minutes and 9 seconds
......@@ -4,6 +4,9 @@
- Small interface change in the template parameters: drop `blocksize` and replace it by `BitSetVector`
- The method can now handle generic `dune-functions` basis types, as long as we have consistency in the data types
- assembleGlobalBasisTransferMatrix:
- Support for all bases in `dune-functions` compatible form added
## Deprecations and removals
- ...
......@@ -7,15 +7,36 @@
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/raviartthomasbasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh>
#include "common.hh"
#include "dune/fufem/test/common.hh"
using namespace Dune;
template<class Matrix, class Function, class FineBasis, class FineCoeff, class CoarseBasis, class CoarseCoeff>
double checkInterpolation(Matrix& matrix, const Function& f,
const FineBasis& fineBasis, FineCoeff& fineCoeffs,
const CoarseBasis& coarseBasis, CoarseCoeff& coarseCoeffs)
{
// assemble the matrix
assembleGlobalBasisTransferMatrix(matrix,coarseBasis,fineBasis);
// interpolate f on the coarse grid
Functions::interpolate(coarseBasis,coarseCoeffs,f);
// interpolate f on the fine grid
Functions::interpolate(fineBasis,fineCoeffs,f);
// now we should have
// fineCoeffs - matrix * coarseCoeffs == 0
matrix.mmv(coarseCoeffs,fineCoeffs);
return fineCoeffs.two_norm2();
}
template<int order, class GW0, class GW1>
bool testLagrange(const GW0& gw0, const GW1& gw1)
{
......@@ -39,22 +60,79 @@ bool testLagrange(const GW0& gw0, const GW1& gw1)
BCRSMatrix<FieldMatrix<double,1,1>> matrix;
BlockVector<FieldVector<double,1>> coarseCoeffs, fineCoeffs;
// assemble the matrix
assembleGlobalBasisTransferMatrix(matrix,coarseLagrange,fineLagrange);
const auto diff = checkInterpolation(matrix, f, fineLagrange, fineCoeffs, coarseLagrange, coarseCoeffs);
// interpolate f on the coarse grid
Functions::interpolate(coarseLagrange,coarseCoeffs,f);
std::cout << " test Lagrange " << GW0::dimension << "D with order " << order << " : error = " << diff << std::endl;
// interpolate f on the fine grid
Functions::interpolate(fineLagrange,fineCoeffs,f);
return diff < 1e-12;
}
// now we should have
// fineCoeffs - matrix * coarseCoeffs == 0
matrix.mmv(coarseCoeffs,fineCoeffs);
std::cout << " test Lagrange " << GW0::dimension << "D with order " << order << " : error = " << fineCoeffs.two_norm2()<< std::endl;
template<int order, class GW0, class GW1>
bool testPowerBasis(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 FieldVector<double,3>{1.,1.,1.};
return fineCoeffs.two_norm2() < 1e-12;
// for the other cases affine function
FieldVector<double,3> y;
y = 1.;
for( auto xs : x )
{
y[0] += xs;
y[1] += 2*xs;
y[2] += 3*xs;
}
return y;
};
using namespace Functions::BasisFactory;
bool passed = true;
{
auto coarseLagrange = makeBasis(gw0, power<3>(lagrange<order>(),blockedInterleaved()));
auto fineLagrange = makeBasis(gw1, power<3>(lagrange<order>(),blockedInterleaved()));
BCRSMatrix<FieldMatrix<double,3,3>> matrix;
BlockVector<FieldVector<double,3>> coarseCoeffs, fineCoeffs;
const auto diff = checkInterpolation(matrix, f, fineLagrange, fineCoeffs, coarseLagrange, coarseCoeffs);
std::cout << " test PowerBasis blockedInterleaved " << GW0::dimension << "D with order " << order << " : error = " << diff << std::endl;
passed = passed and diff < 1e-12;
}
{
auto coarseLagrange = makeBasis(gw0, power<3>(lagrange<order>(),flatInterleaved()));
auto fineLagrange = makeBasis(gw1, power<3>(lagrange<order>(),flatInterleaved()));
BCRSMatrix<FieldMatrix<double,1,1>> matrix;
BlockVector<FieldVector<double,1>> coarseCoeffs, fineCoeffs;
const auto diff = checkInterpolation(matrix, f, fineLagrange, fineCoeffs, coarseLagrange, coarseCoeffs);
std::cout << " test PowerBasis flatInterleaved " << GW0::dimension << "D with order " << order << " : error = " << diff << std::endl;
passed = passed and diff < 1e-12;
}
{
auto coarseLagrange = makeBasis(gw0, power<3>(lagrange<order>(),flatLexicographic()));
auto fineLagrange = makeBasis(gw1, power<3>(lagrange<order>(),flatLexicographic()));
BCRSMatrix<FieldMatrix<double,1,1>> matrix;
BlockVector<FieldVector<double,1>> coarseCoeffs, fineCoeffs;
const auto diff = checkInterpolation(matrix, f, fineLagrange, fineCoeffs, coarseLagrange, coarseCoeffs);
std::cout << " test PowerBasis flatLexicographic " << GW0::dimension << "D with order " << order << " : error = " << diff << std::endl;
passed = passed and diff < 1e-12;
}
return passed;
}
......@@ -156,12 +234,20 @@ bool checkTransferMatrix(const GW0& gw0, const GW1& gw1)
bool passed = true;
// test Lagrange elements for orders 0-3
// for scalar bases ...
passed = passed and testLagrange<0>(gw0,gw1);
passed = passed and testLagrange<1>(gw0,gw1);
passed = passed and testLagrange<2>(gw0,gw1);
// ... and powerBases
passed = passed and testPowerBasis<0>(gw0,gw1);
passed = passed and testPowerBasis<1>(gw0,gw1);
passed = passed and testPowerBasis<2>(gw0,gw1);
// order 3 only for 2D grids
if ( GW0::dimension < 3 )
if constexpr ( GW0::dimension < 3 )
{
passed = passed and testLagrange<3>(gw0,gw1);
passed = passed and testPowerBasis<3>(gw0,gw1);
}
// test transfer operators for Raviart-Thomas bases
if constexpr ( GW0::dimension==2 )
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment