gradientloadmodel.hh

#ifndef __DUNE_ACFEM_MODELS_MODULES_GRADIENTLOADMODEL_HH__
#define __DUNE_ACFEM_MODELS_MODULES_GRADIENTLOADMODEL_HH__
 
#include <dune/fem/function/localfunction/const.hh>
 
#include "../modelbase.hh"
#include "../expressiontraits.hh"
 
namespace Dune {
 
  namespace ACFem::PDEModel {
 
    template<class GridFunction>
    class GradientLoadModel
      : public ModelBase<typename Fem::ToNewDimRangeFunctionSpace<typename std::decay_t<GridFunction>::FunctionSpaceType,
                                                                  std::decay_t<GridFunction>::FunctionSpaceType::dimDomain>::Type>
      , public Expressions::SelfExpression<GradientLoadModel<GridFunction> >
      , public MPL::UniqueTags<ConditionalType<IsConstantExprArg<GridFunction>::value, ConstantExpression, void>,
                               ConditionalType<IsTypedValue<GridFunction>::value, TypedValueExpression, void> >
    {
      static_assert(HasTag<GridFunction, Fem::HasLocalFunction>::value,
                    "GridFunction must provide a local function");
 
      using ThisType = GradientLoadModel;
      using BaseType = ModelBase<typename Fem::ToNewDimRangeFunctionSpace<typename std::decay_t<GridFunction>::FunctionSpaceType,
                                                                          std::decay_t<GridFunction>::FunctionSpaceType::dimDomain>::Type>;
     protected:
      using LocalFunctionType = Fem::ConstLocalFunction<std::decay_t<GridFunction> >;
     public:
      using GridFunctionType = GridFunction;
      using typename BaseType::DomainType;
      using typename BaseType::RangeRangeType;
      using typename BaseType::RangeJacobianRangeType;
 
      using BaseType::dimDomain;
      using BaseType::dimRange;
 
      enum {
        dimWorld = dimDomain
      };
 
      static_assert(dimDomain == dimRange && dimDomain == dimWorld,
                    "This is meant for dimDomain == dimRange, i.e. pressure gradients and such.");
      static_assert(std::decay_t<GridFunction>::FunctionSpaceType::dimRange == 1,
                    "This is only meant for real gradients (not Jacobians), i.e. dimRange of function == 1");
 
      template<class FunctionArg, std::enable_if_t<std::is_constructible<LocalFunctionType, FunctionArg>::value, int> = 0>
      GradientLoadModel(FunctionArg&& function, const std::string& name = "")
        : localFunction_(function)
        , name_(name == "" ? "(grad " + localFunction_.gridFunction().name() + ")" : name)
      {}
 
      std::string name() const
      {
        return name_;
      }
 
      // Interface methods that need to be reimplemented
 
      template<class Entity>
      void bind(const Entity& entity)
      {
        localFunction_.bind(entity);
      }
 
      void unbind()
      {
        localFunction_.unbind();
      }
 
      template<class Quadrature>
      auto flux(const QuadraturePoint<Quadrature> &x) const
      {
        auto scalarValue = localFunction_.evaluate(x);
 
        RangeJacobianRangeType flux(0);
 
        // value * IdentityMatrix
        for (int i = 0; i < dimDomain; ++i) {
          flux[i][i] -= scalarValue;
        }
        return flux;
      }
 
      template<class Quadrature>
      auto fluxDivergence(const QuadraturePoint<Quadrature> &x) const
      {
        return localFunction_.jacobian(x)[0];
      }
 
      template<class Quadrature>
      auto robinFlux(const QuadraturePoint<Quadrature> &x,
                     const DomainType& unitOuterNormal) const
      {
        return RangeRangeType(unitOuterNormal) *= localFunction_.evaluate(x);
      }
 
      template<class Intersection>
      auto classifyBoundary(const Intersection& intersection) const
      {
        // true is correct as higher-level code has to take care of
        // the Dirichlet-(non-Dirichlet) splitting of the boundary.
        return std::make_pair(true, std::bitset<dimRange>());
      }
 
     protected:
      LocalFunctionType localFunction_;
      std::string name_;
    };
 
 
 
    template<class GridFunction>
    constexpr auto
    gradientLoadModel(GridFunction&& f, const std::string& name = "")
    {
      return expressionClosure(GradientLoadModel<GridFunction>(std::forward<GridFunction>(f), name));
    }
 
    template<class GridFunction, std::enable_if_t<ExpressionTraits<GridFunction>::isZero, int> = 0>
    constexpr auto
    gradientLoadModel(GridFunction&& f, const std::string& name = "")
    {
      using FunctionSpace =
        typename Fem::ToNewDimRangeFunctionSpace<typename std::decay_t<GridFunction>::FunctionSpaceType,
                                                 std::decay_t<GridFunction>::FunctionSpaceType::dimDomain>::Type;
 
      return zeroModel(FunctionSpace{}, name);
    }
 
 
 
 
  } // namespace ACFem::PDEModel
 
  namespace ACFem {
 
    using PDEModel::gradientLoadModel;
 
  }
 
}  //Namespace Dune
 
#endif // __DUNE_ACFEM_MODELS_MODULES_GRADIENTLOADMODEL_HH__