Commit c9b8ac08 authored by Patrick Jaap's avatar Patrick Jaap

Merge branch 'feature/functions-mass-matrix-assembler' into 'master'

Implement local mass matrix assembler for arbitrary dune-functions bases

See merge request !79
parents a9b8972a 0724dca9
Pipeline #30847 passed with stage
in 11 minutes and 34 seconds
#pragma once
#include <vector>
#include <dune/common/fvector.hh>
#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
#include <dune/typetree/visitor.hh>
#include <dune/typetree/traversal.hh>
namespace Dune::Fufem
namespace Impl
/** \brief Computes the inner product of the local basis functions
* It computes at each leaf node the inner product (hence, mass matrix) of the
* local basis function. It determines the local index and updates the
* corresponding entry of the local matrix.
template <class LV, class M>
class LeafNodeMassAssembler
using LocalView = LV;
using Matrix = M;
using field_type = typename M::field_type;
/** \brief Constructor with all necessary context information
* \param [in] lv (localView) element information and localBasis
* \param [out] m (matrix) resulting local mass matrix, wrapped in the ISTLBackend
LeafNodeMassAssembler(const LV& lv, M& m)
: localView_(lv)
, matrix_(m)
template<class Node, class TreePath>
void operator()(Node& node, TreePath treePath)
const auto& element = localView_.element();
const auto& geometry = element.geometry();
const auto& finiteElement = node.finiteElement();
const auto& localBasis = finiteElement.localBasis();
constexpr int gridDim = LocalView::GridView::dimension;
using ctype = typename LocalView::GridView::ctype;
std::vector<Dune::FieldVector<field_type,1>> values(localBasis.size());
QuadratureRuleKey quadKey(finiteElement);
const auto& quad = QuadratureRuleCache<double, gridDim>::rule(quadKey);
for( const auto& qp : quad )
const auto& quadPos = qp.position();
const auto integrationElement = geometry.integrationElement(quadPos);
localBasis.evaluateFunction(quadPos, values);
// compute matrix entries = inner products
auto z = qp.weight() * integrationElement;
for(int i=0; i<localBasis.size(); ++i)
double 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)
double zij = values[j] * zi;
auto colIndex = node.localIndex(j);
// mass matrix is symmetric
matrix_[rowIndex][colIndex] += zij;
matrix_[colIndex][rowIndex] += zij;
// z * values[i]^2
matrix_[rowIndex][rowIndex] += values[i] * zi;
const LocalView& localView_;
Matrix& matrix_;
} // namespace Impl
/** \brief Local mass assembler for dune-functions basis
* We assume the mass matrix to be symmetric and hence ansatz and trial space to be equal!
* The implementation allows an arbitrary dune-functions compatible basis, as long as there
* is an ISTLBackend available for the resulting matrix.
class DuneFunctionsLocalMassAssembler
template<class Element, class LocalMatrix, class LocalView>
void operator()(const Element& element, LocalMatrix& localMatrix, const LocalView& trialLocalView, const LocalView& ansatzLocalView) const
// the mass matrix operator assumes trial == ansatz
if (&trialLocalView != &ansatzLocalView)
DUNE_THROW(Dune::Exception,"The mass matrix operator assumes equal ansatz and trial functions");
// matrix was already resized before
localMatrix = 0.;
// create a tree visitor and compute the inner products of the local functions at the leaf nodes
Impl::LeafNodeMassAssembler leadNodeMassAssembler(trialLocalView,localMatrix);
} // namespace Dune::Fufem
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