Commit 2496e2a7 authored by oliver.sander_at_tu-dresden.de's avatar oliver.sander_at_tu-dresden.de

Merge branch 'basistransfermatrix-for-raviart-thomas-elements' into 'master'

Implement and test assembleGlobalBasisTransferMatrix for Raviart-Thomas elements

See merge request !63
parents 191fa340 8d5ca586
Pipeline #27869 failed with stage
in 56 minutes
......@@ -359,9 +359,10 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
// get the local basis
const auto& localBasis = node.finiteElement().localBasis();
using Range = typename std::decay_t<decltype(localBasis)>::Traits::RangeType;
// Hack: The RangeFieldType is the best guess of a suitable type for coefficients we have here
using CoefficientType = typename std::decay_t<decltype(localBasis)>::Traits::RangeFieldType;
std::vector<Range> values(localBasis.size());
std::vector<CoefficientType> coefficients(localBasis.size());
LocalBasisWrapper basisFctj(node.finiteElement().localBasis(),0);
......@@ -376,7 +377,7 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
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);
node.finiteElement().localInterpolation().interpolate(coarseBaseFctj, coefficients);
// get the matrix col index
const auto& globalCoarseIndex = coarseLocalView.index(j);
......@@ -384,7 +385,8 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
// set the matrix indices
for (size_t i=0; i<localBasis.size(); i++) {
if ( std::abs(values[i]) < tolerance )
// some values may be practically zero -- no need to store those
if ( std::abs(coefficients[i]) < tolerance )
continue;
// get the matrix row index
......@@ -446,9 +448,9 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
// get the local basis
const auto& localBasis = node.finiteElement().localBasis();
using Range = typename std::decay_t<decltype(localBasis)>::Traits::RangeType;
using Coefficient = typename std::decay_t<decltype(localBasis)>::Traits::RangeFieldType;
std::vector<Range> values(localBasis.size());
std::vector<Coefficient> coefficients(localBasis.size());
LocalBasisWrapper basisFctj(node.finiteElement().localBasis(),0);
......@@ -463,7 +465,7 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
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);
node.finiteElement().localInterpolation().interpolate(coarseBaseFctj, coefficients);
// get the matrix col index
const auto& globalCoarseIndex = coarseLocalView.index(j);
......@@ -471,7 +473,8 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
// set the matrix indices
for (size_t i=0; i<localBasis.size(); i++) {
if ( std::abs(values[i]) < tolerance )
// some values may be practically zero -- no need to store those
if ( std::abs(coefficients[i]) < tolerance )
continue;
// get the matrix row index
......@@ -480,7 +483,7 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
if (processed[globalFineIndex][0])
continue;
Dune::MatrixVector::addToDiagonal(matrix[globalFineIndex][globalCoarseIndex],values[i]);
Dune::MatrixVector::addToDiagonal(matrix[globalFineIndex][globalCoarseIndex],coefficients[i]);
}
}
......
......@@ -6,6 +6,7 @@
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/raviartthomasbasis.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
......@@ -51,22 +52,77 @@ bool testLagrange(const GW0& gw0, const GW1& gw1)
// fineCoeffs - matrix * coarseCoeffs == 0
matrix.mmv(coarseCoeffs,fineCoeffs);
std::cout << " test lagrage " << GW0::dimension << "D with order " << order << " : error = " << fineCoeffs.two_norm2()<< std::endl;
std::cout << " test Lagrange " << 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>
bool testRaviartThomas(const GW0& gw0, const GW1& gw1)
{
std::cout << " test Raviart-Thomas bases " << GW0::dimension << "D" << std::endl;
Functions::RaviartThomasBasis<GW0,order> coarseBasis(gw0);
Functions::RaviartThomasBasis<GW1,order> fineBasis(gw1);
// TODO Remove FieldMatrix/FieldVector here!
BCRSMatrix<FieldMatrix<double,1,1>> matrix;
BlockVector<FieldVector<double,1>> coarseCoeffs, fineCoeffs;
// assemble the matrix
assembleGlobalBasisTransferMatrix(matrix,coarseBasis,fineBasis);
// construct some nontrivial function in the coarse-grid Raviart-Thomas space
coarseCoeffs.resize(coarseBasis.dimension());
for (std::size_t i=0; i<coarseCoeffs.size(); i++)
coarseCoeffs[i] = i*std::sqrt(2.0);
// use the transfer matrix to inject that function into the fine-grid space
fineCoeffs.resize(fineBasis.dimension());
matrix.mv(coarseCoeffs,fineCoeffs);
// the two space are nested, hence both functions should be the same.
// Let's check that!
using Range = FieldVector<double,GW0::dimension>;
auto coarseFunction = Functions::makeDiscreteGlobalBasisFunction<Range>(coarseBasis, coarseCoeffs);
auto fineFunction = Functions::makeDiscreteGlobalBasisFunction<Range>( fineBasis, fineCoeffs);
auto coarseLocalFunction = localFunction(coarseFunction);
auto fineLocalFunction = localFunction(fineFunction);
for (const auto& element : elements(gw1))
{
coarseLocalFunction.bind(element.father());
fineLocalFunction.bind(element);
FieldVector<double,2> testPoint(0);
auto testPointInFather = element.geometryInFather().global(testPoint);
// TODO: The run-time test is disabled for the time being, because it requires
// https://gitlab.dune-project.org/staging/dune-functions/-/merge_requests/241
if ( false && (coarseLocalFunction(testPointInFather) - fineLocalFunction(testPoint) ).two_norm() > 1e-10 )
{
std::cerr << "Raviart-Thomas: transfer matrix is not the natural injection!" << std::endl;
return false;
}
}
return true;
}
template<int order, class GW0, class GW1, class GW2>
bool testTwoLevelLagrange(const GW0& gw0, const GW1& gw1, const GW2& gw2)
{
// create 3 global basis
// create 3 global bases
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
// matrices for the combined jump and two simple jumps
BCRSMatrix<FieldMatrix<double,1,1>> matrixTwoJumps, matrixJumpOne, matrixJumpTwo;
// assemble the matrices
......@@ -74,7 +130,7 @@ bool testTwoLevelLagrange(const GW0& gw0, const GW1& gw1, const GW2& gw2)
assembleGlobalBasisTransferMatrix(matrixJumpOne,coarseLagrange,semiLagrange);
assembleGlobalBasisTransferMatrix(matrixJumpTwo,semiLagrange,fineLagrange);
// now check wether the product of two jumps is the conbined jump
// now check whether the product of two jumps is the combined jump
// compute matrixTwoJumps -= matrixJumpTwo*matrixJumpOne
for (auto row0 = matrixJumpTwo.begin(); row0 != matrixJumpTwo.end(); row0++)
for (auto entry0 = row0->begin(); entry0 != row0->end(); entry0++)
......@@ -93,13 +149,13 @@ bool testTwoLevelLagrange(const GW0& gw0, const GW1& gw1, const GW2& gw2)
// this tests the interpolation porperty of the matrix
// this tests the interpolation property of the matrix
template<class GW0, class GW1>
bool checkTransferMatrix(const GW0& gw0, const GW1& gw1)
{
bool passed = true;
// test oder 0-3
// test Lagrange elements for orders 0-3
passed = passed and testLagrange<0>(gw0,gw1);
passed = passed and testLagrange<1>(gw0,gw1);
passed = passed and testLagrange<2>(gw0,gw1);
......@@ -107,6 +163,10 @@ bool checkTransferMatrix(const GW0& gw0, const GW1& gw1)
if ( GW0::dimension < 3 )
passed = passed and testLagrange<3>(gw0,gw1);
// test transfer operators for Raviart-Thomas bases
if constexpr ( GW0::dimension==2 )
passed = passed and testRaviartThomas<0>(gw0,gw1);
return passed;
}
......
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