basisinterpolationmatrixassembler.hh

Go to the documentation of this file.

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set ts=8 sw=2 et sts=2:
 
// 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_BASIS_INTERPOLATION_MATRIX_ASSEMBLER_HH
#define DUNE_FUFEM_ASSEMBLERS_BASIS_INTERPOLATION_MATRIX_ASSEMBLER_HH
 
#include <type_traits>
#include <cmath>
#include <vector>
 
#include <dune/common/exceptions.hh>
#include <dune/common/version.hh>
#include <dune/common/rangeutilities.hh>
 
#include <dune/assembler/backends/istlmatrixbackend.hh>
 
#if DUNE_VERSION_GTE(DUNE_FUNCTIONS, 2, 11)
#include <dune/common/typetree/childaccess.hh>
#include <dune/common/typetree/traversal.hh>
#else
#include <dune/typetree/childextraction.hh>
#include <dune/typetree/traversal.hh>
#endif
 
 
 
namespace Impl{
 
 
// This class provides the local interpolation operator
// from a coarse local finite element to a fine local finite
// element. The associated grid elements are mapped using a geometry
// map. After being created, this can be bound to the coarse
// basis functions and afterwards provides the respective interpolation
// coefficients.
//
// Notice that this is restricted to cases, where FineLocalInterpolation::interpolate()
// evaluates at a fixed deterministic point set not depending on the passed function.
// In particular this will not work if interpolate() uses adaptive quadrature
// or if it evaluatuated derivatives.
template<class CoarseFiniteElement, class FineFiniteElement, class Geometry>
class LocalInterpolationOperator
{
  using CoarseLocalBasis = typename CoarseFiniteElement::Traits::LocalBasisType;
  using CoarseRange = typename CoarseLocalBasis::Traits::RangeType;
 
  // Hack: The RangeFieldType is the best guess of a suitable type for coefficients we have here
  using CoefficientType = typename CoarseLocalBasis::Traits::RangeFieldType;
 
  using FineLocalInterpolation = typename FineFiniteElement::Traits::LocalInterpolationType;
 
public:
 
  LocalInterpolationOperator(const CoarseFiniteElement& coarseFiniteElement, const FineFiniteElement& fineFiniteElement, const Geometry& geometry)
    : coarseLocalBasis_(coarseFiniteElement.localBasis())
    , fineLocalInterpolation_(fineFiniteElement.localInterpolation())
    , geometry_(geometry)
    , coefficients_(fineFiniteElement.size())
    , coarseValues_(fineFiniteElement.size())
  {
    // We first compute the values at all interpolation points.
    // To this end a dummy call to interpolate() is used to deduce
    // the interpolation points and evaluate the coarse bases.
    // Later the, we can simply return the pre-computed values
    // based on the evaluation order.
    // This requires that interpolate() does the evaluation at
    // a deterministic point set in a deterministic order and
    // thus rules out e.g. adaptive quadrature.
    // If the interpolation is given by simple point evaluations,
    // the cache computed here is essentially the local interpolation
    // matrix. However, it is also fine, if interpolate() does
    // additional computations based on the values, like, e.g.
    // a Piola-transformation.
    auto k = 0;
    auto trackingCoarseFunction = [&](const auto& x) {
      coarseValues_.resize(k+1);
      coarseLocalBasis_.evaluateFunction(geometry_(x), coarseValues_[k]);
      return coarseValues_[k++][0];
    };
    fineLocalInterpolation_.interpolate(trackingCoarseFunction, coefficients_);
  }
 
  const auto& bindToCoarseBasisFunction(std::size_t coarseIndex)
  {
    // baseFctj behaves like the j-th coarse bases function if
    // interpolate() evaluates it at the same points in the same
    // order as above.
    auto k = 0;
    auto baseFctj = [&](const auto& x) {
      return coarseValues_[k++][coarseIndex];
    };
 
    // Interpolate j^th base function by the fine basis
    fineLocalInterpolation_.interpolate(baseFctj, coefficients_);
    return coefficients_;
  }
 
  const auto& coefficients() const
  {
    return coefficients_;
  }
 
private:
  const CoarseLocalBasis& coarseLocalBasis_;
  const FineLocalInterpolation& fineLocalInterpolation_;
  const Geometry& geometry_;
  std::vector<CoefficientType> coefficients_;
  std::vector<std::vector<CoarseRange>> coarseValues_;
};
 
 
 
} // namespace Impl
 
 
template<class MatrixType, class CoarseBasis, class FineBasis>
static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
                                              const CoarseBasis& coarseBasis,
                                              const FineBasis& fineBasis,
                                              const double tolerance = 1e-12)
{
  const auto& coarseGridView = coarseBasis.gridView();
  const auto& fineGridView   = fineBasis.gridView();
 
  auto coarseLocalView = coarseBasis.localView();
  auto fineLocalView   = fineBasis.localView();
 
  // for multi index access wrap the matrix to the ISTLBackend
  constexpr bool isAssemblerMatrixBackend = requires(MatrixType& matrixBackend) {
    MatrixType::LocalPattern(0,0);
    matrixBackend.setZero();
  };
 
  constexpr bool isFufemMatrixBackend = requires(MatrixType& matrixBackend) {
    matrixBackend.matrix();
    matrixBackend.assign(0);
  };
 
  // For compatibility we support three modes:
  // Passing a plain matrix, a dune-fufem backend,
  // or a dune-assembler backend.
  auto&& matrixBackend = [&]() -> decltype(auto) {
    if constexpr(isAssemblerMatrixBackend)
      return matrix;
    else if constexpr(isFufemMatrixBackend)
      return Dune::Assembler::ISTLMatrixBackend(matrix.matrix());
    else
      return Dune::Assembler::ISTLMatrixBackend(matrix);
  }();
 
  // prepare index set with multi index access
  auto patternBuilder = matrixBackend.patternBuilder();
  patternBuilder.resize(fineBasis, coarseBasis);
 
  // ///////////////////////////////////////////
  // Determine the indices present in the matrix
  // ///////////////////////////////////////////
 
  // This map transforms coordinates to the ancestor element
  // by chaining the GeometryInFather entries in the vector.
  using GeometryInFather = std::decay_t<decltype(elements(fineGridView).begin()->geometryInFather())>;
  std::vector<GeometryInFather> geometryInFathersVector;
  const auto geometryInFathers = [&](auto&& x) {
    auto y = x;
    for (const auto& g : geometryInFathersVector)
      y = g.global(y);
    return y;
  };
 
  using MatrixBackend = std::decay_t<decltype(matrixBackend)>;
 
  auto localPattern = typename MatrixBackend::LocalPattern(fineLocalView.maxSize(), coarseLocalView.maxSize());
 
  // loop over all fine grid elements
  for ( const auto& fineElement : elements(fineGridView) )
  {
    fineLocalView.bind(fineElement);
 
    // find a coarse element that is the parent of the fine element
    // start in the fine grid and track the geometry mapping
    auto coarseElement = fineElement;
    // collect all geometryInFather's on the way
    geometryInFathersVector.clear();
    while ( not coarseGridView.contains(coarseElement) and coarseElement.hasFather() )
    {
      // add the geometry to the container
      geometryInFathersVector.push_back(coarseElement.geometryInFather());
      // step up one level
      coarseElement = coarseElement.father();
    }
 
    // did we really end up in coarseElement?
    if ( not coarseGridView.contains(coarseElement) )
      DUNE_THROW( Dune::Exception, "There is a fine element without a parent in the coarse grid!");
 
    coarseLocalView.bind(coarseElement);
 
    localPattern.resize(fineLocalView.size(), coarseLocalView.size());
    localPattern.clear();
 
    Dune::TypeTree::forEachLeafNode(coarseLocalView.tree(), [&] (const auto& coarseNode, const auto& treePath) {
      auto&& fineNode = Dune::TypeTree::child(fineLocalView.tree(), treePath);
 
      auto localInterpolationOperator = Impl::LocalInterpolationOperator(coarseNode.finiteElement(), fineNode.finiteElement(), geometryInFathers);
      const auto& coefficients = localInterpolationOperator.coefficients();
 
      for (auto j : Dune::range(coarseNode.size()))
      {
        localInterpolationOperator.bindToCoarseBasisFunction(j);
        for (auto i : Dune::range(fineNode.size()))
        {
          if (std::abs(coefficients[i]) >= tolerance)
            localPattern.add(fineNode.localIndex(i), coarseNode.localIndex(j));
        }
      }
    });
    patternBuilder.scatter(fineLocalView, coarseLocalView, localPattern);
  }
 
  // Setup matrix from computed index set
  patternBuilder.setupMatrix();
 
  // Now set the values
 
  for ( const auto& fineElement : elements(fineGridView) )
  {
    fineLocalView.bind(fineElement);
 
    // find a coarse element that is the parent of the fine element
    // start in the fine grid and track the geometry mapping
    auto coarseElement = fineElement;
    // collect all geometryInFather's on the way
    geometryInFathersVector.clear();
    while ( not coarseGridView.contains(coarseElement) and coarseElement.hasFather() )
    {
      // add the geometry to the container
      geometryInFathersVector.push_back(coarseElement.geometryInFather());
      // step up one level
      coarseElement = coarseElement.father();
    }
 
    coarseLocalView.bind(coarseElement);
 
    // Ideally we would assemble into a local matrix and use scatter,
    // but the latter only allows to add values. Here we want to
    // overwrite them.
    Dune::TypeTree::forEachLeafNode(coarseLocalView.tree(), [&] (const auto& coarseNode, const auto& treePath) {
      auto&& fineNode = Dune::TypeTree::child(fineLocalView.tree(), treePath);
 
      auto localInterpolationOperator = Impl::LocalInterpolationOperator(coarseNode.finiteElement(), fineNode.finiteElement(), geometryInFathers);
      const auto& coefficients = localInterpolationOperator.coefficients();
 
      for (auto j : Dune::range(coarseNode.size()))
      {
        localInterpolationOperator.bindToCoarseBasisFunction(j);
        const auto& globalCoarseIndex = coarseLocalView.index( coarseNode.localIndex( j ) );
        for (auto i : Dune::range(fineNode.size()))
        {
          if (std::abs(coefficients[i]) >= tolerance)
          {
            const auto& globalFineIndex = fineLocalView.index( fineNode.localIndex( i ) );
            matrixBackend(globalFineIndex, globalCoarseIndex) = coefficients[i];
          }
        }
      }
    });
 
  }
}
 
 
#endif