weakdivergenceloadmodel.hh

#ifndef __DUNE_ACFEM_MODELS_MODULES_WEAKDIVERGENCELOADMODEL_HH__
#define __DUNE_ACFEM_MODELS_MODULES_WEAKDIVERGENCELOADMODEL_HH__
 
#include <dune/fem/function/localfunction/const.hh>
 
#include "../../expressions/terminal.hh"
#include "../../functions/functions.hh"
#include "../modelbase.hh"
 
namespace Dune {
 
  namespace ACFem::PDEModel {
 
    template<class GridFunction>
    class WeakDivergenceLoadModel
      : public ModelBase<typename std::decay_t<GridFunction>::FunctionSpaceType::ScalarFunctionSpaceType>
      , public Expressions::SelfExpression<WeakDivergenceLoadModel<GridFunction> >
      , public MPL::UniqueTags<ConditionalType<ExpressionTraits<GridFunction>::isVolatile, VolatileExpression, void>,
                               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 = WeakDivergenceLoadModel;
      using BaseType = ModelBase<typename std::decay_t<GridFunction>::FunctionSpaceType::ScalarFunctionSpaceType>;
      using GridFunctionDecay = std::decay_t<GridFunction>;
      using LocalFunctionType = Fem::ConstLocalFunction<GridFunctionDecay>;
     public:
      using GridFunctionType = GridFunction;
 
      using typename BaseType::RangeType;
      using typename BaseType::DomainType;
      using typename BaseType::JacobianRangeType;
      using BaseType::dimDomain;
      using BaseType::dimRange;
 
      enum {
        dimWorld = GridFunctionDecay::GridPartType::dimensionworld
      };
 
      static_assert(dimRange == 1 && dimDomain == dimWorld,
                    "This is only meant for dimRange == 1 and dimDomain == dimWorld.");
 
      template<class FunctionArg, std::enable_if_t<std::is_constructible<LocalFunctionType, FunctionArg>::value, int> = 0>
      WeakDivergenceLoadModel(FunctionArg&& function, const std::string& name = "")
        : localFunction_(std::forward<FunctionArg>(function))
        , name_(name == "" ? "(div " + 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>
      JacobianRangeType flux(const QuadraturePoint<Quadrature> &x) const
      {
        const DomainType dataValue = localFunction_.evaluate(x);
 
        JacobianRangeType flux = 0;
        // value * IdentityMatrix
        for (int i = 0; i < dimRange; ++i) {
          flux[i] -= dataValue;
        }
        return flux;
      }
 
      template<class Quadrature>
      auto fluxDivergence(const QuadraturePoint<Quadrature> &x) const
      {
        const auto divergenceData = localFunction_.jacobian(x);
 
        RangeType result = 0;
        for (int i = 0; i < dimWorld; ++i) {
          result[0] += divergenceData[i][i];
        }
        return result;
      }
 
     protected:
      LocalFunctionType localFunction_;
      std::string name_;
    };
 
 
 
    template<class GridFunction>
    constexpr auto weakDivergenceModel(GridFunction&& f, const std::string& name = "")
    {
      static_assert(IsWrappableByConstLocalFunction<GridFunction>::value,
                    "GridFunction is not suitable for Fem::ConstLocalFunction");
 
      return expressionClosure(WeakDivergenceLoadModel<GridFunction>(std::forward<GridFunction>(f), name));
    }
 
    template<class GridFunction, std::enable_if_t<ExpressionTraits<GridFunction>::isZero, int> = 0>
    constexpr auto weakDivergenceModel(GridFunction&& f, const std::string& name = "")
    {
      using FunctionSpaceType = typename std::decay_t<GridFunction>::FunctionSpaceType::ScalarFunctionSpaceType;
      return zeroModel(FunctionSpaceType{});
    }
 
 
 
 
  } // namespace ACFem
 
}  //Namespace Dune
 
 
#endif // __DUNE_ACFEM_MODELS_MODULES_WEAKDIVERGENCELOADMODEL_HH__