Fix and test assembleGlobalBasisTransferMatrix for Raviart-Thomas elements

The run-time test in this patch is disabled, because it requires

  https://gitlab.dune-project.org/staging/dune-functions/-/merge_requests/241

which is not merged yet.  However, Patrick Jaap is currently
improving the assembleGlobalBasisTransferMatrix method,
and I want the Raviart-Thomas compile-time check in as quickly
as possible, so he can you it to test his changes.
parent de8dfb3f
Pipeline #27868 failed with stage
in 63 minutes and 47 seconds
......@@ -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>
......@@ -57,6 +58,61 @@ bool testLagrange(const GW0& gw0, const GW1& gw1)
}
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)
{
......@@ -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