Skip to content
Snippets Groups Projects

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

Merged Patrick Jaap requested to merge fix/blocked-local-matrix into master
All threads resolved!
1 file
+ 20
9
Compare changes
  • Side-by-side
  • Inline
@@ -62,6 +62,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 +144,19 @@ 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::FieldMatrix<double,dim,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!
Dune::Matrix<Dune::FieldMatrix<double,dim,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
for(size_t i=0; i<localMatrix.N(); i++)
for(size_t j=0; j<localMatrix.M(); j++)
localMatrix[i][j] = localBlockedMatrix[i/dim][j/dim][i%dim][j%dim];
};
@@ -185,13 +192,19 @@ 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::FieldMatrix<double,dim,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!
Dune::Matrix<Dune::FieldMatrix<double,dim,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
for(size_t i=0; i<localMatrix.N(); i++)
for(size_t j=0; j<localMatrix.M(); j++)
localMatrix[i][j] = localBlockedMatrix[i/dim][j/dim][i%dim][j%dim];
};
auto massMatrixBackend = Dune::Fufem::istlMatrixBackend(localMassMatrix);
@@ -256,7 +269,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 +294,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(),
Loading