Skip to content
Snippets Groups Projects
Commit 8ef53aa6 authored by oliver.sander_at_tu-dresden.de's avatar oliver.sander_at_tu-dresden.de
Browse files

Fix mass assembler for vector-valued bases

And add a unit test.
parent 1d36a4c8
No related branches found
No related tags found
1 merge request!90Make mass assembler work with vector-valued bases
......@@ -48,10 +48,10 @@ public:
const auto& geometry = element.geometry();
const auto& finiteElement = node.finiteElement();
const auto& localBasis = finiteElement.localBasis();
using RangeType = typename std::decay_t<decltype(localBasis)>::Traits::RangeType;
constexpr int gridDim = LocalView::GridView::dimension;
using ctype = typename LocalView::GridView::ctype;
std::vector<Dune::FieldVector<field_type,1>> values(localBasis.size());
std::vector<RangeType> values(localBasis.size());
QuadratureRuleKey quadKey(finiteElement);
const auto& quad = QuadratureRuleCache<double, gridDim>::rule(quadKey);
......@@ -65,16 +65,16 @@ public:
// compute matrix entries = inner products
auto z = qp.weight() * integrationElement;
for(int i=0; i<localBasis.size(); ++i)
for (std::size_t i=0; i<localBasis.size(); ++i)
{
double zi = values[i]*z;
auto zi = values[i]*z;
auto rowIndex = node.localIndex(i);
// start off-diagonal, treat the diagonal extra
for (int j=i+1; j<localBasis.size(); ++j)
for (std::size_t j=i+1; j<localBasis.size(); ++j)
{
double zij = values[j] * zi;
auto zij = values[j] * zi;
auto colIndex = node.localIndex(j);
......
......@@ -22,6 +22,7 @@ set(GRID_BASED_TESTS
integraloperatorassemblertest
istlbackendtest
laplaceassemblertest
massassemblertest
polynomialtest
secondorderassemblertest
subgridxyfunctionalassemblertest
......
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set ts=8 sw=2 et sts=2:
#include <config.h>
#include <array>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/raviartthomasbasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/fufem/assemblers/dunefunctionsoperatorassembler.hh>
#include <dune/fufem/assemblers/istlbackend.hh>
#include <dune/fufem/assemblers/localassemblers/dunefunctionsmassassembler.hh>
using namespace Dune;
template <class Basis, class Matrix>
void test(const Basis& basis)
{
using Assembler = Fufem::DuneFunctionsOperatorAssembler<Basis, Basis>;
auto matrix = Matrix{};
auto matrixBackend = Dune::Fufem::istlMatrixBackend(matrix);
auto assembler = Assembler{basis, basis};
auto patternBuilder = matrixBackend.patternBuilder();
assembler.assembleBulkPattern(patternBuilder);
patternBuilder.setupMatrix();
Fufem::DuneFunctionsLocalMassAssembler localMassAssembler;
assembler.assembleBulkEntries(matrixBackend, localMassAssembler);
}
int main (int argc, char *argv[])
{
Dune::MPIHelper::instance(argc, argv);
const int dim = 2;
// Build a test grid
using Grid = YaspGrid<dim>;
Grid grid({1.0, 1.0}, {5, 5});
auto gridView = grid.leafGridView();
using namespace Dune::Functions::BasisBuilder;
// Test with a scalar bases
{
auto basis = makeBasis(
gridView,
lagrange<2>()
);
using Basis = decltype(basis);
using Matrix = BCRSMatrix<double>;
test<Basis,Matrix>(basis);
}
// Test with a basis of vector-valued functions
{
auto basis = makeBasis(
gridView,
raviartThomas<1>()
);
using Basis = decltype(basis);
using Matrix = BCRSMatrix<double>;
test<Basis,Matrix>(basis);
}
// Test with a power basis
{
auto basis = makeBasis(
gridView,
power<2>(
lagrange<2>(),
flatInterleaved()
));
using Basis = decltype(basis);
using Matrix = BCRSMatrix<double>;
test<Basis,Matrix>(basis);
}
return 0;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment