globalvaluedlocalfiniteelement.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
// SPDX-FileCopyrightText: Copyright © DUNE 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_FUNCTIONS_FUNCTIONSPACEBASES_GLOBALVALUEDLOCALFINITEELEMENT_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_GLOBALVALUEDLOCALFINITEELEMENT_HH
 
#include <array>
#include <numeric>
 
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/math.hh>
#include <dune/common/rangeutilities.hh>
 
#include <dune/geometry/referenceelements.hh>
 
#include <dune/localfunctions/common/localbasis.hh>
#include <dune/localfunctions/common/localfiniteelementtraits.hh>
 
#include <dune/functions/common/densevectorview.hh>
 
namespace Dune::Functions::Impl
{
 
  struct ContravariantPiolaTransformator
  {
    template<typename Values, typename LocalCoordinate, typename Geometry>
    static auto apply(Values& values,
                      const LocalCoordinate& xi,
                      const Geometry& geometry)
    {
      auto jacobianTransposed = geometry.jacobianTransposed(xi);
      auto integrationElement = geometry.integrationElement(xi);
 
      for (auto& value : values)
      {
        auto tmp = value;
        // Here we cannot directly use
        // jacobianTransposed.mtv(tmp, value);
        // because mtv expects the DenseVector interface.
        auto tmpDenseVector = Impl::DenseVectorView(tmp);
        auto valueDenseVector = Impl::DenseVectorView(value);
        jacobianTransposed.mtv(tmpDenseVector, valueDenseVector);
        value /= integrationElement;
      }
    }
 
    template<typename Gradients, typename LocalCoordinate, typename Geometry>
    static auto applyJacobian(Gradients& gradients,
                              const LocalCoordinate& xi,
                              const Geometry& geometry)
    {
      auto jacobianTransposed = geometry.jacobianTransposed(xi);
      auto integrationElement = geometry.integrationElement(xi);
      for (auto& gradient : gradients)
      {
        auto tmp = gradient;
        gradient = 0;
        for (size_t k=0; k<gradient.M(); k++)
          for (size_t l=0; l<tmp.N(); l++)
            // Use sparseRange because jacobianTransposed may be a sparse DiagonalMatrix
            for(auto&& [jacobianTransposed_l_j, j] : sparseRange(jacobianTransposed[l]))
              gradient[j][k] += jacobianTransposed_l_j * tmp[l][k];
        gradient /= integrationElement;
      }
    }
 
    template<class Function, class LocalCoordinate, class Element>
    class LocalValuedFunction
    {
      const Function& f_;
      const Element& element_;
 
      using LocalValue = LocalCoordinate;
 
    public:
 
      LocalValuedFunction(const Function& f, const Element& element)
      : f_(f), element_(element)
      {}
 
      auto operator()(const LocalCoordinate& xi) const
      {
        auto globalValue = f_(xi);
 
        // Apply the inverse Piola transform
        auto jacobianInverseTransposed = element_.geometry().jacobianInverseTransposed(xi);
        auto integrationElement = element_.geometry().integrationElement(xi);
 
        auto localValue = LocalValue{};
 
        // Here we cannot directly use
        // jacobianInverseTransposed.mtv(globalValue, localValue);
        // because mtv expects the DenseVector interface.
        auto globalValueDenseVector = Impl::DenseVectorView(globalValue);
        auto localValueDenseVector = Impl::DenseVectorView(localValue);
        jacobianInverseTransposed.mtv(globalValueDenseVector, localValueDenseVector);
        localValue *= integrationElement;
 
        return localValue;
      }
    };
  };
 
  struct CovariantPiolaTransformator
  {
    template<typename Values, typename LocalCoordinate, typename Geometry>
    static auto apply(Values& values,
                      const LocalCoordinate& xi,
                      const Geometry& geometry)
    {
      auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(xi);
 
      for (auto& value : values)
      {
        auto tmp = value;
        // Here we cannot directly use
        // jacobianInverseTransposed.mv(tmp, value);
        // because mv expects the DenseVector interface.
        auto tmpDenseVector = Impl::DenseVectorView(tmp);
        auto valueDenseVector = Impl::DenseVectorView(value);
        jacobianInverseTransposed.mv(tmpDenseVector, valueDenseVector);
      }
    }
 
    template<typename Gradients, typename LocalCoordinate, typename Geometry>
    static auto applyJacobian(Gradients& gradients,
                              const LocalCoordinate& xi,
                              const Geometry& geometry)
    {
      auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(xi);
 
      for (auto& gradient : gradients)
      {
        auto tmp = gradient;
        gradient = 0;
        for (size_t j=0; j<gradient.N(); j++)
          for (size_t k=0; k<gradient.M(); k++)
            // Use sparseRange because jacobianTransposed may be a sparse DiagonalMatrix
            for(auto&& [jacobianInverseTransposed_j_l, l] : sparseRange(jacobianInverseTransposed[j]))
              gradient[j][k] += jacobianInverseTransposed_j_l * tmp[l][k];
      }
    }
 
    template<class Function, class LocalCoordinate, class Element>
    class LocalValuedFunction
    {
      const Function& f_;
      const Element& element_;
 
    public:
 
      LocalValuedFunction(const Function& f, const Element& element)
      : f_(f), element_(element)
      {}
 
      auto operator()(const LocalCoordinate& xi) const
      {
        auto globalValue = f_(xi);
 
        // Apply the inverse Piola transform
        auto jacobianTransposed = element_.geometry().jacobianTransposed(xi);
 
        auto localValue = globalValue;
        // Here we cannot directly use
        // jacobianTransposed.mv(globalValue, localValue);
        // because mv expects the DenseVector interface.
        auto localValueDenseVector = Impl::DenseVectorView(localValue);
        auto globalValueDenseVector = Impl::DenseVectorView(globalValue);
        jacobianTransposed.mv(globalValueDenseVector, localValueDenseVector);
 
        return localValue;
      }
    };
  };
 
  template<class Transformator, class LocalValuedLocalBasis, class Element>
  class GlobalValuedLocalBasis
  {
  public:
    using Traits = typename LocalValuedLocalBasis::Traits;
 
    void bind(const LocalValuedLocalBasis& localValuedLocalBasis, const Element& element)
    {
      localValuedLocalBasis_ = &localValuedLocalBasis;
      element_ = &element;
    }
 
    auto size() const
    {
      return localValuedLocalBasis_->size();
    }
 
    void evaluateFunction(const typename Traits::DomainType& x,
                          std::vector<typename Traits::RangeType>& out) const
    {
      localValuedLocalBasis_->evaluateFunction(x,out);
 
      Transformator::apply(out, x, element_->geometry());
    }
 
    void evaluateJacobian(const typename Traits::DomainType& x,
                          std::vector<typename Traits::JacobianType>& out) const
    {
      localValuedLocalBasis_->evaluateJacobian(x,out);
 
      Transformator::applyJacobian(out, x, element_->geometry());
    }
 
    void partial(const std::array<unsigned int,2>& order,
                 const typename Traits::DomainType& x,
                 std::vector<typename Traits::RangeType>& out) const
    {
      auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
      if (totalOrder == 0) {
        evaluateFunction(x, out);
      } else if (totalOrder == 1) {
        auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
        out.resize(size());
 
        // TODO: The following is wasteful:  We compute the full Jacobian and then return
        // only a part of it.  While we need the full Jacobian of the underlying local-valued LFE,
        // it should be possible to compute only a partial Piola transform for the requested
        // partial derivatives.
        std::vector<typename Traits::JacobianType> fullJacobian;
        localValuedLocalBasis_->evaluateJacobian(x,fullJacobian);
 
        Transformator::applyJacobian(fullJacobian, x, element_->geometry());
 
        for (std::size_t i=0; i<out.size(); i++)
          for (std::size_t j=0; j<out[i].size(); j++)
            out[i][j] = fullJacobian[i][j][direction];
 
      } else
        DUNE_THROW(NotImplemented, "Partial derivatives of order 2 or higher");
    }
 
    auto order() const
    {
      return localValuedLocalBasis_->order();
    }
 
    const LocalValuedLocalBasis* localValuedLocalBasis_;
    const Element* element_;
  };
 
  template<class Transformator, class LocalValuedLocalInterpolation, class Element>
  class GlobalValuedLocalInterpolation
  {
  public:
    void bind(const LocalValuedLocalInterpolation& localValuedLocalInterpolation, const Element& element)
    {
      localValuedLocalInterpolation_ = &localValuedLocalInterpolation;
      element_ = &element;
    }
 
    template<typename F, typename C>
    void interpolate (const F& f, std::vector<C>& out) const
    {
      using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
      typename Transformator::template LocalValuedFunction<F,LocalCoordinate,Element> localValuedFunction(f, *element_);
      localValuedLocalInterpolation_->interpolate(localValuedFunction, out);
    }
 
  private:
    const LocalValuedLocalInterpolation* localValuedLocalInterpolation_;
    const Element* element_;
  };
 
 
  template<class Transformator, class LocalValuedLFE, class Element>
  class GlobalValuedLocalFiniteElement
  {
    using LocalBasis = GlobalValuedLocalBasis<Transformator,
                                              typename LocalValuedLFE::Traits::LocalBasisType,
                                              Element>;
    using LocalInterpolation = GlobalValuedLocalInterpolation<Transformator,
                                                              typename LocalValuedLFE::Traits::LocalInterpolationType,
                                                              Element>;
 
  public:
    using Traits = LocalFiniteElementTraits<LocalBasis,
                                            typename LocalValuedLFE::Traits::LocalCoefficientsType,
                                            LocalInterpolation>;
 
    GlobalValuedLocalFiniteElement() {}
 
    void bind(const LocalValuedLFE& localValuedLFE, const Element& element)
    {
      globalValuedLocalBasis_.bind(localValuedLFE.localBasis(), element);
      globalValuedLocalInterpolation_.bind(localValuedLFE.localInterpolation(), element);
      localValuedLFE_ = &localValuedLFE;
    }
 
    const typename Traits::LocalBasisType& localBasis() const
    {
      return globalValuedLocalBasis_;
    }
 
    const typename Traits::LocalCoefficientsType& localCoefficients() const
    {
      return localValuedLFE_->localCoefficients();
    }
 
    const typename Traits::LocalInterpolationType& localInterpolation() const
    {
      return globalValuedLocalInterpolation_;
    }
 
    std::size_t size() const
    {
      return localValuedLFE_->size();
    }
 
    GeometryType type() const
    {
      return localValuedLFE_->type();
    }
 
  private:
 
    typename Traits::LocalBasisType globalValuedLocalBasis_;
    typename Traits::LocalInterpolationType globalValuedLocalInterpolation_;
    const LocalValuedLFE* localValuedLFE_;
  };
 
}        // namespace Dune::Functions::Impl
 
#endif   // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_GLOBALVALUEDLOCALFINITEELEMENT_HH