diff --git a/dune/fufem/assemblers/localassemblers/dunefunctionsmassassembler.hh b/dune/fufem/assemblers/localassemblers/dunefunctionsmassassembler.hh new file mode 100644 index 0000000000000000000000000000000000000000..9868e76d78be2e95835cd8dfd8da9ee8219ee0cb --- /dev/null +++ b/dune/fufem/assemblers/localassemblers/dunefunctionsmassassembler.hh @@ -0,0 +1,126 @@ +#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 +{ +public: + + 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; + } + } + } + +private: + 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 +{ +public: + + 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); + Dune::TypeTree::forEachLeafNode(trialLocalView.tree(),leadNodeMassAssembler); + } + +}; + +} // namespace Dune::Fufem