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

Fix: TrustRegionSolver: use blocked matrices for fufem scalar basis assembler

valgrind complained here and it is right: in the internal loop of the
assembler we loop over localMatrix.N() and M(). Since dune-functions
uses flat matrices for the localMatrix the number of dof's is dim-times
too large for a powerBasis. In order to spare a double implementation of
the local assembler only for powerBasis, we can simply convert it in the
lambda. valgrind is now happy here.
parent 46ec7d63
No related branches found
No related tags found
1 merge request!31Fix: TrustRegionSolver: use blocked matrices for fufem scalar basis assembler
Pipeline #28509 passed
......@@ -2,6 +2,7 @@
#include <dune/common/timer.hh>
#include <dune/istl/io.hh>
#include <dune/istl/scaledidmatrix.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
......@@ -62,6 +63,7 @@ setup(const typename BasisType::GridView::Grid& grid,
damping_ = damping;
int numLevels = grid_->maxLevel()+1;
const auto dim = VectorType::value_type::dimension;
#if HAVE_DUNE_PARMG
Dune::Timer setupTimer;
......@@ -143,13 +145,21 @@ setup(const typename BasisType::GridView::Grid& grid,
{
using FiniteElement = std::decay_t<decltype(basis.localView().tree().child(0).finiteElement())>;
// construct a Fufem Basis Assembler
auto laplaceStiffness = LaplaceAssembler<GridType, FiniteElement, FiniteElement>();
auto laplaceStiffness = LaplaceAssembler<GridType, FiniteElement, FiniteElement, Dune::ScaledIdentityMatrix<double,dim>>();
// transform it to a dune-functions assembler
auto localAssembler = [&](const auto& element, auto& localMatrix, auto&& trialLocalView, auto&& ansatzLocalView){
laplaceStiffness.assemble(element,
localMatrix,
trialLocalView.tree().child(0).finiteElement(),
ansatzLocalView.tree().child(0).finiteElement());
// the fufem basis assembler assumes a scalar basis but a blocked element matrix!
// we can use ScaledIdentityMatrix here
Dune::Matrix<Dune::ScaledIdentityMatrix<double,dim>> localBlockedMatrix( localMatrix.N()/dim, localMatrix.M()/dim );
laplaceStiffness.assemble(element,
localBlockedMatrix,
trialLocalView.tree().child(0).finiteElement(),
ansatzLocalView.tree().child(0).finiteElement());
// convert back to flat matrix (we need only the diagonal entries in the blocks)
localMatrix = 0;
for(size_t i=0; i<localBlockedMatrix.N(); i++)
for(size_t j=0; j<localBlockedMatrix.M(); j++)
Dune::MatrixVector::addToDiagonal(localMatrix[i*dim][j*dim], localBlockedMatrix[i][j].scalar());
};
......@@ -185,13 +195,21 @@ setup(const typename BasisType::GridView::Grid& grid,
{
// construct a Fufem Basis Assembler
using FiniteElement = std::decay_t<decltype(basis.localView().tree().child(0).finiteElement())>;
auto massStiffness = MassAssembler<GridType, FiniteElement, FiniteElement>();
auto massStiffness = MassAssembler<GridType, FiniteElement, FiniteElement, Dune::ScaledIdentityMatrix<double,dim>>();
// transform it to a dune-functions assembler
auto localMassAssembler = [&](const auto& element, auto& localMatrix, auto&& trialLocalView, auto&& ansatzLocalView){
// the fufem basis assembler assumes a scalar basis but a blocked element matrix!
// we can use ScaledIdentityMatrix here
Dune::Matrix<Dune::ScaledIdentityMatrix<double,dim>> localBlockedMatrix( localMatrix.N()/dim, localMatrix.M()/dim );
massStiffness.assemble(element,
localMatrix,
localBlockedMatrix,
trialLocalView.tree().child(0).finiteElement(),
ansatzLocalView.tree().child(0).finiteElement());
// convert back to flat matrix (we need only the diagonal entries in the blocks)
localMatrix = 0;
for(size_t i=0; i<localBlockedMatrix.N(); i++)
for(size_t j=0; j<localBlockedMatrix.M(); j++)
Dune::MatrixVector::addToDiagonal(localMatrix[i*dim][j*dim], localBlockedMatrix[i][j].scalar());
};
auto massMatrixBackend = Dune::Fufem::istlMatrixBackend(localMassMatrix);
......@@ -256,7 +274,6 @@ setup(const typename BasisType::GridView::Grid& grid,
if constexpr ( Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() )
{
const auto dim = VectorType::value_type::dimension;
using namespace Dune::Functions::BasisFactory;
auto p1PowerBasis = makeBasis(
basis.gridView(),
......@@ -282,7 +299,6 @@ setup(const typename BasisType::GridView::Grid& grid,
if constexpr ( Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() )
{
const auto dim = VectorType::value_type::dimension;
using namespace Dune::Functions::BasisFactory;
auto p1PowerBasis = makeBasis(
basis.gridView(),
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment