laplaceassembler.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_LAPLACEASSEMBLER_HH
#define DUNE_FUFEM_ASSEMBLERS_LOCALASSEMBLERS_LAPLACEASSEMBLER_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 <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/version.hh>
 
#include <dune/istl/matrix.hh>
 
#if DUNE_VERSION_GTE(DUNE_FUNCTIONS, 2, 11)
#include <dune/common/typetree/traversal.hh>
#else
#include <dune/typetree/traversal.hh>
#endif
 
#include "dune/fufem/quadraturerules/quadraturerulecache.hh"
 
 
 
namespace Dune::Fufem
{
 
  namespace Impl
  {
 
    template <class LV, class M>
    class LeafNodeLaplaceAssembler
    {
    public:
 
      using LocalView = LV;
      using Matrix = M;
      using field_type = typename M::field_type;
 
      LeafNodeLaplaceAssembler(const LV& lv, M& m)
      : localView_(lv)
      , matrix_(m)
      {}
 
      template<class Node, class TreePath>
      void operator()(Node& node, [[maybe_unused]] TreePath treePath)
      {
        const auto& element = localView_.element();
        const auto& geometry = element.geometry();
        const auto& finiteElement = node.finiteElement();
        const auto& localBasis = finiteElement.localBasis();
        using JacobianType = typename std::decay_t<decltype(localBasis)>::Traits::JacobianType;
        constexpr int gridDim = LocalView::GridView::dimension;
 
        auto feQuadKey = QuadratureRuleKey(finiteElement);
        // we have a product of jacobians
        auto quadKey = feQuadKey.derivative().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);
 
          std::vector<JacobianType> referenceJacobians;
          localBasis.evaluateJacobian(quadPos, referenceJacobians);
 
          const auto& jacobianInverse = geometry.jacobianInverse(quadPos);
 
          // transform jacobians
          std::vector<decltype(referenceJacobians[0] * jacobianInverse)> jacobians(referenceJacobians.size());
          for (size_t i=0; i<jacobians.size(); i++)
            jacobians[i] = referenceJacobians[i] * jacobianInverse;
 
          // compute matrix entries = inner products of jacobians
          auto z = qp.weight() * integrationElement;
          for (std::size_t i=0; i<localBasis.size(); ++i)
            for (std::size_t j=0; j<localBasis.size(); ++j)
              for ( std::size_t ii=0; ii<jacobians[i].N(); ii++ )
                for ( std::size_t jj=0; jj<jacobians[i].M(); jj++ )
                  matrix_[node.localIndex(i)][node.localIndex(j)] += z*jacobians[i][ii][jj]*jacobians[j][ii][jj];
        }
      }
 
    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.")]]
  LaplaceAssembler
  {
  public:
 
    [[deprecated("This class is deprecated and will be removed after 2.11. Use Dune::Fufem::Forms instead.")]]
    LaplaceAssembler()
    {}
 
    template<class Element, class LocalMatrix, class LocalView>
    void operator()(const Element& element, LocalMatrix& localMatrix, const LocalView& testLocalView, const LocalView& ansatzLocalView) const
    {
      // the Laplace matrix operator assumes test == ansatz
      if (&testLocalView != &ansatzLocalView)
        DUNE_THROW(Dune::Exception,"The Laplace 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::LeafNodeLaplaceAssembler leadNodeLaplaceAssembler(testLocalView,localMatrix);
      TypeTree::forEachLeafNode(testLocalView.tree(),leadNodeLaplaceAssembler);
    }
 
  };
 
} // namespace Dune::Fufem
 
 
 
#endif  // DUNE_FUFEM_ASSEMBLERS_LOCALASSEMBLERS_LAPLACEASSEMBLER_HH