massassembler.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUFEM_ASSEMBLERS_LOCALASSEMBLERS_MASSASSEMBLER_HH
#define DUNE_FUFEM_ASSEMBLERS_LOCALASSEMBLERS_MASSASSEMBLER_HH
 
#ifndef DUNE_FUFEM_DISABLE_HEADER_DEPRECATION
#warning This header is deprecated and will be removed after 2.11. Use Dune::Fufem::Forms instead.
#endif
 
#include <vector>
 
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/typetree/traversal.hh>
 
#include <dune/istl/matrix.hh>
 
#include "dune/fufem/quadraturerules/quadraturerulecache.hh"
 
 
namespace Dune::Fufem
{
 
  namespace Impl
  {
 
    template <class LV, class M>
    class LeafNodeMassAssembler
    {
    public:
 
      using LocalView = LV;
      using Matrix = M;
      using field_type = typename M::field_type;
 
      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();
        using RangeType = typename std::decay_t<decltype(localBasis)>::Traits::RangeType;
        constexpr int gridDim = LocalView::GridView::dimension;
 
        std::vector<RangeType> values(localBasis.size());
 
        auto feQuadKey = QuadratureRuleKey(finiteElement);
        // we have a product of FE functions
        auto quadKey = feQuadKey.square();
        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 (std::size_t i=0; i<localBasis.size(); ++i)
          {
            auto zi = values[i]*z;
 
            auto rowIndex = node.localIndex(i);
 
            // start off-diagonal, treat the diagonal extra
            for (std::size_t j=i+1; j<localBasis.size(); ++j)
            {
              auto 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
 
  class
  [[deprecated("This class is deprecated and will be removed after 2.11. Use Dune::Fufem::Forms instead.")]]
  MassAssembler
  {
  public:
 
    [[deprecated("This class is deprecated and will be removed after 2.11. Use Dune::Fufem::Forms instead.")]]
    MassAssembler()
    {}
 
    template<class Element, class LocalMatrix, class LocalView>
    void operator()(const Element& element, LocalMatrix& localMatrix, const LocalView& testLocalView, const LocalView& ansatzLocalView) const
    {
      // the mass matrix operator assumes test == ansatz
      if (&testLocalView != &ansatzLocalView)
        DUNE_THROW(Dune::Exception,"The mass matrix operator assumes equal ansatz and test functions");
 
      // create a tree visitor and compute the inner products of the local functions at the leaf nodes
      Impl::LeafNodeMassAssembler leadNodeMassAssembler(testLocalView,localMatrix);
      TypeTree::forEachLeafNode(testLocalView.tree(),leadNodeMassAssembler);
    }
 
  };
 
} // namespace Dune::Fufem
 
#endif  // DUNE_FUFEM_ASSEMBLERS_LOCALASSEMBLERS_MASSASSEMBLER_HH