transformedoperator.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 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_FORMS_TRANSFORMEDOPERATOR_HH
#define DUNE_FUFEM_FORMS_TRANSFORMEDOPERATOR_HH
 
#include <cstddef>
#include <type_traits>
#include <tuple>
#include <utility>
#include <functional>
 
#include <dune/fufem/forms/baseclass.hh>
#include <dune/fufem/forms/tensors.hh>
#include <dune/fufem/forms/localfunctionadaptor.hh>
#include <dune/fufem/forms/quadraturehints.hh>
 
 
 
namespace Dune::Fufem::Forms::Impl {
 
  template<class Operator, int arity = Operator::arity>
  struct TransformedOperatorBaseTraits
  {
    using type = MultilinearOperator<arity>;
  };
 
  template<class Operator>
  struct TransformedOperatorBaseTraits<Operator, 1>
  {
    using type = UnaryOperator<Operator::argIndex>;
  };
 
} // namespace Dune::Fufem::Forms::Impl
 
 
 
namespace Dune::Fufem::Forms {
 
 
 
  template<class Op, class BaseOperator, class QuadratureHint=Dune::Fufem::Forms::QuadratureHints::KeepOrder>
  class TransformedOperator
    : public Impl::TransformedOperatorBaseTraits<BaseOperator>::type
  {
    using Transformation = Op;
    using BaseLocalOperator = decltype(localOperator(std::declval<BaseOperator>()));
 
  public:
 
    using Element = typename BaseOperator::Element;
    using Range = decltype(std::declval<Op>()(std::declval<typename BaseOperator::Range>()));
 
    TransformedOperator(const Transformation& transformation, const BaseOperator& baseOperator, QuadratureHint quadratureHint)
      : baseOperator_(baseOperator)
      , transformation_(transformation)
      , quadratureHint_(quadratureHint)
    {}
 
    TransformedOperator(const Transformation& transformation, const BaseOperator& baseOperator)
    requires std::is_default_constructible_v<QuadratureHint>
      : baseOperator_(baseOperator)
      , transformation_(transformation)
      , quadratureHint_()
    {}
 
    TransformedOperator(const BaseOperator& baseOperator)
    requires std::is_default_constructible_v<Transformation> && std::is_default_constructible_v<QuadratureHint>
      : baseOperator_(baseOperator)
      , transformation_()
      , quadratureHint_()
    {}
 
 
    class LocalOperator : public BaseLocalOperator
    {
    public:
 
      using Element = typename BaseLocalOperator::Element;
      using Intersection = typename BaseLocalOperator::Intersection;
      using Range = typename TransformedOperator::Range;
 
      LocalOperator(const Transformation& transformation, BaseLocalOperator&& baseLocalOperator, QuadratureHint quadratureHint)
        : BaseLocalOperator(std::move(baseLocalOperator))
        , transformation_(transformation)
        , quadratureHint_(quadratureHint)
      {}
 
      void bind(const Element& element)
      {
        BaseLocalOperator::bind(element);
        quadratureRuleKey_ = quadratureHint_(BaseLocalOperator::quadratureRuleKey());
      }
 
      void bind(const Intersection& intersection, const Element& element, const Element& otherElement)
      {
        BaseLocalOperator::bind(intersection, element, otherElement);
        quadratureRuleKey_ = quadratureHint_(BaseLocalOperator::quadratureRuleKey());
      }
 
      auto quadratureRuleKey() const
      {
        return quadratureRuleKey_;
      }
 
      auto operator()(std::size_t index) const
      {
        return compose(std::cref(transformation_), BaseLocalOperator::operator()(index));
      }
 
    private:
      Transformation transformation_;
      QuadratureHint quadratureHint_;
      QuadratureRuleKey quadratureRuleKey_;
    };
 
    friend LocalOperator localOperator(const TransformedOperator& transformedOperator)
    {
      return LocalOperator(
        transformedOperator.transformation(),
        localOperator(transformedOperator.baseOperator()),
        transformedOperator.quadratureHint_
      );
    }
 
    friend LocalFunctionAdaptor<LocalOperator> localFunction(const TransformedOperator& transformedOperator)
    requires (TransformedOperator::arity==0)
    {
      return LocalFunctionAdaptor<LocalOperator>(localOperator(transformedOperator));
    }
 
    auto basis() const
    {
      return baseOperator_.basis();
    }
 
    auto treePath() const
    {
      return baseOperator_.treePath();
    }
 
    const BaseOperator& baseOperator() const
    {
      return baseOperator_;
    }
 
    const Transformation& transformation() const
    {
      return transformation_;
    }
 
    friend auto gradient(const TransformedOperator& f)
    {
      static constexpr auto dimension = Element::Geometry::Coordinate::dimension;
      static_assert(dimension<4);
 
      auto gradientTransformation = [transformation = f.transformation()](const auto& g) {
        static_assert(dimension<=3, "The gradient of a transformed operator is not implemented for range dimensions > 3.");
        if constexpr (dimension==1)
          return std::array{transformation(g[0])};
        else if constexpr (dimension==2)
          return std::array{transformation(g[0]), transformation(g[1])};
        else if constexpr (dimension==3)
          return std::array{transformation(g[0]), transformation(g[1]), transformation(g[2])};
      };
      return compose(gradientTransformation, gradient(f.baseOperator()));
    }
 
    friend auto grad(const TransformedOperator& f)
    {
      return gradient(f);
    }
 
    template<std::size_t k>
    auto rebindArgIndex() const
    {
      auto reboundBaseOperator = baseOperator_.template rebindArgIndex<k>();
      using ReboundBaseOperator = decltype(reboundBaseOperator);
      return TransformedOperator<Op, ReboundBaseOperator>(transformation_, reboundBaseOperator);
    }
 
  private:
    BaseOperator baseOperator_;
    Transformation transformation_;
    QuadratureHint quadratureHint_;
  };
 
 
 
} // namespace Dune::Fufem::Forms
 
 
#endif // DUNE_FUFEM_FORMS_TRANSFORMEDOPERATOR_HH