combinedoperator.hh

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=8 sw=2 sts=2:
#ifndef DUNE_PDELAB_LOCALOPERATOR_COMBINEDOPERATOR_HH
#define DUNE_PDELAB_LOCALOPERATOR_COMBINEDOPERATOR_HH
 
#include <cstddef>
 
#include <tuple>
#include <dune/common/tupleutility.hh>
#include <dune/common/typetraits.hh>
 
#include <dune/pdelab/localoperator/callswitch.hh>
 
namespace Dune {
  namespace PDELab {
 
 
 
    template<typename ApplyOp, typename... Args>
    class CombinedOperator
    {
    protected:
      using ArgPtrs = std::tuple<std::shared_ptr<std::remove_reference_t<Args>>...>;
      using ArgRefs = std::tuple<Args&...>;
 
      ArgPtrs lops;
 
      template<typename... FArgs>
      void applyLops(FArgs &... args) const
      {
        static_cast<const ApplyOp&>(*this).applyLops(args...);
      }
 
    public:
      //
      //
 
      CombinedOperator() {}
 
      CombinedOperator(Args&... args) : lops(stackobject_to_shared_ptr(args)...) {}
 
      CombinedOperator(Args&&... args) : lops(std::make_shared<Args>(std::move(args))...) {}
 
    protected:
      // CombinedOperator(const ArgPtrs & l) : lops(l) {}
      CombinedOperator(ArgPtrs && l) : lops(std::move(l)) {}
    public:
 
      template<std::size_t i>
      void setSummand(typename std::tuple_element_t<i,ArgRefs> summand)
      { std::get<i>(lops) = &summand; }
 
      template<std::size_t i>
      typename std::tuple_element_t<i,ArgRefs> getSummand()
      { return *std::get<i>(lops); }
 
 
      //
      //
 
    private:
      // we only evaluate the one-/two-sided flags, if a skeleton term is to be evaluated
      template<typename T>
      using OneSidedSkeleton = std::integral_constant
      < bool, ( ( T::doAlphaSkeleton || T::doLambdaSkeleton) && ! T::doSkeletonTwoSided)>;
      template<typename T>
      using TwoSidedSkeleton = std::integral_constant
      < bool, ( ( T::doAlphaSkeleton || T::doLambdaSkeleton) && T::doSkeletonTwoSided)>;
 
    public:
 
      enum { doSkipEntity            =
             std::disjunction_v<std::integral_constant<bool,Args::doSkipEntity>...>                   };
 
 
      enum { doSkipIntersection      =
             std::disjunction_v<std::integral_constant<bool,Args::doSkipIntersection>...>             };
 
 
      enum { doPatternVolume             =
             std::disjunction_v<std::integral_constant<bool,Args::doPatternVolume>...>                };
 
      enum { doPatternVolumePostSkeleton =
             std::disjunction_v<std::integral_constant<bool,Args::doPatternVolumePostSkeleton>...>    };
      enum { doPatternSkeleton           =
             std::disjunction_v<std::integral_constant<bool,Args::doPatternSkeleton>...>              };
      enum { doPatternBoundary           =
             std::disjunction_v<std::integral_constant<bool,Args::doPatternBoundary>...>              };
 
 
      enum { doAlphaVolume               =
             std::disjunction_v<std::integral_constant<bool,Args::doAlphaVolume>...>                  };
      enum { doAlphaVolumePostSkeleton   =
             std::disjunction_v<std::integral_constant<bool,Args::doAlphaVolumePostSkeleton>...>      };
      enum { doAlphaSkeleton             =
             std::disjunction_v<std::integral_constant<bool,Args::doAlphaSkeleton>...>                };
      enum { doAlphaBoundary             =
             std::disjunction_v<std::integral_constant<bool,Args::doAlphaBoundary>...>                };
 
      enum { doLambdaVolume              =
             std::disjunction_v<std::integral_constant<bool,Args::doLambdaVolume>...>                 };
      enum { doLambdaVolumePostSkeleton  =
             std::disjunction_v<std::integral_constant<bool,Args::doLambdaVolumePostSkeleton>...>     };
      enum { doLambdaSkeleton            =
             std::disjunction_v<std::integral_constant<bool,Args::doLambdaSkeleton>...>               };
      enum { doLambdaBoundary            =
             std::disjunction_v<std::integral_constant<bool,Args::doLambdaBoundary>...>               };
 
      enum { doSkeletonTwoSided          =
             std::disjunction_v<TwoSidedSkeleton<Args>...>     };
      static_assert(!(std::conjunction_v<OneSidedSkeleton<Args>...> &&
                      std::conjunction_v<TwoSidedSkeleton<Args>...>),
                      "Some summands require a one-sided skelton, others a "
                      "two-sided skeleton.  This is not supported.");
 
      enum { isLinear = std::conjunction_v<std::integral_constant<bool,Args::isLinear>...> };
 
 
      //
      //
 
      template<typename EG>
      bool skip_entity
      ( const EG& eg) const
      {
        // check that all operators return the same
        bool different_skips = false;
        Hybrid::forEach(std::make_index_sequence<sizeof...(Args)>{},
          [&](auto i) {
            different_skips ^= getSummand<i>().skip_entity(eg);
            }
          );
        if (different_skips)
          DUNE_THROW(RangeError, "CombinedOperator is not allowed to have "
                                  "different skip_entity results");
        return getSummand<0>().skip_entity(eg);
      }
 
      template<typename IG>
      bool skip_intersection
      ( const IG& ig) const
      {
        // check that all operators return the same
        bool different_skips = false;
        Hybrid::forEach(std::make_index_sequence<sizeof...(Args)>{},
          [&](auto i) {
            different_skips ^= getSummand<i>().skip_intersection(ig);
            }
          );
        if (different_skips)
          DUNE_THROW(RangeError, "CombinedOperator is not allowed to have "
                                  "different skip_intersection results");
        return getSummand<0>().skip_intersection(ig);
      }
 
 
      //
      //
 
 
      template<typename LFSU, typename LFSV, typename LocalPattern>
      void pattern_volume
      ( const LFSU& lfsu, const LFSV& lfsv,
        LocalPattern& pattern) const
      {
        applyLops(LocalOperatorApply::patternVolume, lfsu, lfsv, pattern);
      }
 
 
      template<typename LFSU, typename LFSV, typename LocalPattern>
      void pattern_volume_post_skeleton
      ( const LFSU& lfsu, const LFSV& lfsv,
        LocalPattern& pattern) const
      {
        applyLops(LocalOperatorApply::patternVolumePostSkeleton, lfsu, lfsv, pattern);
      }
 
 
      template<typename LFSU, typename LFSV, typename LocalPattern>
      void pattern_skeleton
      ( const LFSU& lfsu_s, const LFSV& lfsv_s,
        const LFSU& lfsu_n, const LFSV& lfsv_n,
        LocalPattern& pattern_sn,
        LocalPattern& pattern_ns) const
      {
        applyLops(LocalOperatorApply::patternSkeleton,
          lfsu_s, lfsv_s, lfsu_n, lfsv_n,
          pattern_sn, pattern_ns);
      }
 
 
      template<typename LFSU, typename LFSV, typename LocalPattern>
      void pattern_boundary
      ( const LFSU& lfsu_s, const LFSV& lfsv_s,
        LocalPattern& pattern_ss) const
      {
        applyLops(LocalOperatorApply::patternBoundary, lfsu_s, lfsv_s, pattern_ss);
      }
 
 
      //
      //
 
 
      template<typename EG, typename LFSU, typename X, typename LFSV,
               typename R>
      void alpha_volume
      ( const EG& eg,
        const LFSU& lfsu, const X& x, const LFSV& lfsv,
        R& r) const
      {
        applyLops(LocalOperatorApply::alphaVolume, eg, lfsu, x, lfsv, r);
      }
 
 
      template<typename EG, typename LFSU, typename X, typename LFSV,
               typename R>
      void alpha_volume_post_skeleton
      ( const EG& eg,
        const LFSU& lfsu, const X& x, const LFSV& lfsv,
        R& r) const
      {
        applyLops(LocalOperatorApply::alphaVolumePostSkeleton, eg, lfsu, x, lfsv, r);
      }
 
 
      template<typename IG, typename LFSU, typename X, typename LFSV,
               typename R>
      void alpha_skeleton
      ( const IG& ig,
        const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
        const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
        R& r_s, R& r_n) const
      {
        applyLops(LocalOperatorApply::alphaSkeleton, ig,
                lfsu_s, x_s, lfsv_s,
                lfsu_n, x_n, lfsv_n,
                r_s, r_n);
      }
 
 
      template<typename IG, typename LFSU, typename X, typename LFSV,
               typename R>
      void alpha_boundary
      ( const IG& ig,
        const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
        R& r_s) const
      {
        applyLops(LocalOperatorApply::alphaVolumePostSkeleton, ig, lfsu_s, x_s, lfsv_s, r_s);
      }
 
 
      //
      //
 
 
      template<typename EG, typename LFSV, typename R>
      void lambda_volume(const EG& eg, const LFSV& lfsv, R& r) const
      {
        applyLops(LocalOperatorApply::lambdaVolume, eg, lfsv, r);
      }
 
 
      template<typename EG, typename LFSV, typename R>
      void lambda_volume_post_skeleton(const EG& eg,
                                       const LFSV& lfsv,
                                       R& r) const
      {
        applyLops(LocalOperatorApply::lambdaVolumePostSkeleton, eg, lfsv, r);
      }
 
 
      template<typename IG, typename LFSV, typename R>
      void lambda_skeleton(const IG& ig,
                           const LFSV& lfsv_s, const LFSV& lfsv_n,
                           R& r_s, R& r_n) const
      {
        applyLops(LocalOperatorApply::lambdaSkeleton, ig, lfsv_s, lfsv_n, r_s, r_n);
      }
 
 
      template<typename IG, typename LFSV, typename R>
      void lambda_boundary(const IG& ig, const LFSV& lfsv_s, R& r_s) const
      {
        applyLops(LocalOperatorApply::lambdaBoundary, ig, lfsv_s, r_s);
      }
 
 
      //
      //
 
 
      template<typename EG, typename LFSU, typename X, typename LFSV,
               typename Y>
      void jacobian_apply_volume
      ( const EG& eg,
        const LFSU& lfsu, const X& x, const LFSV& lfsv,
        Y& y) const
      {
        applyLops(LocalOperatorApply::jacobianApplyVolume, eg, lfsu, x, lfsv, y);
      }
 
 
      template<typename EG, typename LFSU, typename X, typename LFSV,
               typename Y>
      void jacobian_apply_volume_post_skeleton
      ( const EG& eg,
        const LFSU& lfsu, const X& x, const LFSV& lfsv,
        Y& y) const
      {
        applyLops(LocalOperatorApply::jacobianApplyVolumePostSkeleton, eg, lfsu, x, lfsv, y);
      }
 
 
      template<typename IG, typename LFSU, typename X, typename LFSV,
               typename Y>
      void jacobian_apply_skeleton
      ( const IG& ig,
        const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
        const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
        Y& y_s, Y& y_n) const
      {
        applyLops(LocalOperatorApply::jacobianApplySkeleton, ig,
                lfsu_s, x_s, lfsv_s,
                lfsu_n, x_n, lfsv_n,
                y_s, y_n);
      }
 
 
      template<typename IG, typename LFSU, typename X, typename LFSV,
               typename Y>
      void jacobian_apply_boundary
      ( const IG& ig,
        const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
        Y& y_s) const
      {
        applyLops(LocalOperatorApply::jacobianApplyBoundary, ig, lfsu_s, x_s, lfsv_s, y_s);
      }
 
 
      //
      //
 
 
      template<typename EG, typename LFSU, typename X, typename LFSV,
               typename LocalMatrix>
      void jacobian_volume
      ( const EG& eg,
        const LFSU& lfsu, const X& x, const LFSV& lfsv,
        LocalMatrix& mat) const
      {
        applyLops(LocalOperatorApply::jacobianVolume, eg, lfsu, x, lfsv, mat);
      }
 
 
      template<typename EG, typename LFSU, typename X, typename LFSV,
               typename LocalMatrix>
      void jacobian_volume_post_skeleton
      ( const EG& eg,
        const LFSU& lfsu, const X& x, const LFSV& lfsv,
        LocalMatrix& mat) const
      {
        applyLops(LocalOperatorApply::jacobianVolumePostSkeleton, eg, lfsu, x, lfsv, mat);
      }
 
 
      template<typename IG, typename LFSU, typename X, typename LFSV,
               typename LocalMatrix>
      void jacobian_skeleton
      ( const IG& ig,
        const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
        const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
        LocalMatrix& mat_ss, LocalMatrix& mat_sn,
        LocalMatrix& mat_ns, LocalMatrix& mat_nn) const
      {
        applyLops(LocalOperatorApply::jacobianSkeleton, ig,
                lfsu_s, x_s, lfsv_s,
                lfsu_n, x_n, lfsv_n,
                mat_ss, mat_sn, mat_ns, mat_nn);
      }
 
 
      template<typename IG, typename LFSU, typename X, typename LFSV,
               typename LocalMatrix>
      void jacobian_boundary
      ( const IG& ig,
        const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
        LocalMatrix& mat_ss) const
      {
        applyLops(LocalOperatorApply::jacobianBoundary, ig, lfsu_s, x_s, lfsv_s, mat_ss);
      }
 
 
      //
      //
 
      typedef typename std::tuple_element<0, std::tuple<Args...>>::type::RealType RealType;
 
    private:
      // template meta program helpers for the methods related to instationary
      // stuff
 
      struct Apply {
      template<typename LOP>
      static void setTime(const LOP& lop, RealType t)
      {
        lop.setTime(t);
      }
 
      template<typename LOP>
      static void preStep(const LOP& lop,
                          RealType time, RealType dt, int stages)
      {
        lop.preStep(time, dt, stages);
      }
 
      template<typename LOP>
      static void postStep(const LOP& lop)
      {
        lop.postStep();
      }
 
      template<typename LOP>
      static void preStage(const LOP& lop, RealType time, int r)
      {
        lop.preStage(time, r);
      }
 
      template<typename LOP>
      static void postStage(const LOP& lop)
      {
        lop.postStage();
      }
 
      template<typename LOP>
      static RealType suggestTimestep(const LOP& lop, RealType & dt)
      {
        dt = std::min(dt,lop.suggestTimestep(dt));
      }
      };
 
    public:
      void setTime (RealType t)
      {
        applyLops(Apply::setTime, t);
      }
 
      RealType getTime () const
      {
        return get<0>(lops)->getTime();
      }
 
      void preStep (RealType time, RealType dt, int stages)
      {
        applyLops(Apply::preStep, time, dt, stages);
      }
 
      void postStep ()
      {
        applyLops(Apply::postStep, lops);
      }
 
      void preStage (RealType time, int r)
      {
        applyLops(Apply::preStage, time, r);
      }
 
      int getStage () const
      {
        return get<0>(lops)->getStage();
      }
 
      void postStage ()
      {
        applyLops(Apply::postStage, lops);
      }
 
 
      RealType suggestTimestep (RealType dt) const
      {
        applyLops(Apply::suggestTimestep, dt);
        return dt;
      }
 
    };
 
  }
}
 
#endif // DUNE_PDELAB_LOCALOPERATOR_COMBINEDOPERATOR_HH