neumannmodel.hh

#ifndef __DUNE_ACFEM_MODELS_MODULES_NEUMANNMODEL_HH__
#define __DUNE_ACFEM_MODELS_MODULES_NEUMANNMODEL_HH__
 
#include "../../expressions/expressionoperations.hh"
#include "../../functions/localfunctiontraits.hh"
#include "../../functions/functiontraits.hh"
 
#include "../indicators/boundaryindicator.hh"
#include "../modelbase.hh"
 
namespace Dune {
 
  namespace ACFem::PDEModel {
 
    template<class GridFunction, class Indicator = EntireBoundaryIndicator>
    class NeumannBoundaryModel
      : public ModelBase<typename std::decay_t<GridFunction>::FunctionSpaceType>
      , public Expressions::SelfExpression<NeumannBoundaryModel<GridFunction, Indicator> >
      , public MPL::UniqueTags<ConditionalType<ExpressionTraits<GridFunction>::isVolatile
                                               ||
                                               ExpressionTraits<Indicator>::isVolatile,
                                               VolatileExpression, void>,
                               ConditionalType<IsConstantExprArg<GridFunction>::value
                                               &&
                                               IsConstantExprArg<Indicator>::value,
                                               ConstantExpression, void>,
                               ConditionalType<IsTypedValue<GridFunction>::value
                                               &&
                                               IsTypedValue<Indicator>::value,
                                               TypedValueExpression, void> >
    {
      static_assert(HasTag<GridFunction, Fem::HasLocalFunction>::value,
                    "GridFunction must provide a local function");
 
      using ThisType = NeumannBoundaryModel;
      using BaseType = ModelBase<typename std::decay_t<GridFunction>::FunctionSpaceType>;
    protected:
      using LocalFunctionType = Fem::ConstLocalFunction<std::decay_t<GridFunction> >;
     public:
      using GridFunctionType = GridFunction;
      using IndicatorType = Indicator;
      using DomainType = typename BaseType::DomainType;
      using RangeType = typename BaseType::RangeType;
      using BoundaryConditionsType = typename BaseType::BoundaryConditionsType;
      using BaseType::dimRange;
 
      template<class FunctionArg, class IndicatorArg,
               std::enable_if_t<(std::is_constructible<LocalFunctionType, FunctionArg>::value
                                 && std::is_constructible<Indicator, IndicatorArg>::value
                                 && !std::is_default_constructible<IndicatorArg>::value
                 ), int> = 0>
      NeumannBoundaryModel(FunctionArg&& values,
                           IndicatorArg&& indicator,
                           const std::string& name = "")
        : localFunction_(std::forward<FunctionArg>(values))
        , indicator_(std::forward<IndicatorArg>(indicator))
        , supported_(false, std::bitset<dimRange>())
        , name_(name == "" ? "Neumann(" + localFunction_.gridFunction().name() + ")" : name)
      {}
 
      template<class FunctionArg, class IndicatorArg,
               std::enable_if_t<(std::is_constructible<LocalFunctionType, FunctionArg>::value
                                 && std::is_constructible<Indicator, IndicatorArg>::value
                                 && std::is_default_constructible<IndicatorArg>::value
                 ), int> = 0>
      NeumannBoundaryModel(FunctionArg&& values,
                           IndicatorArg&& indicator = IndicatorArg(),
                           const std::string& name = "")
        : localFunction_(std::forward<FunctionArg>(values))
        , indicator_(std::forward<IndicatorArg>(indicator))
        , supported_(false, std::bitset<dimRange>())
        , name_(name == "" ? "Neumann(" + localFunction_.gridFunction().name() + ")" : name)
      {}
 
      std::string name() const
      {
        return name_;
      }
 
      template<class Entity>
      void bind(const Entity& entity)
      {
        if constexpr (!IndicatorTraits<Indicator>::emptySupport) {
          localFunction_.bind(entity);
        }
      }
 
      void unbind()
      {
        if constexpr (!IndicatorTraits<Indicator>::emptySupport) {
          localFunction_.unbind();
        }
      }
 
      template<class Intersection>
      auto classifyBoundary(const Intersection& intersection)
      {
        supported_.first = IndicatorTraits<IndicatorType>::globalSupport || indicator_.applies(intersection);
        return supported_;
      }
 
      template<class Quadrature>
      RangeType robinFlux(const QuadraturePoint<Quadrature> &x,
                          const DomainType& unitOuterNormal) const
      {
        if (!IndicatorTraits<Indicator>::emptySupport
            &&
            (IndicatorTraits<Indicator>::globalSupport || supported_.first)) {
          RangeType result;
          localFunction_.evaluate(x, result);
          return -result;
        }
        return RangeType(0.);
      }
 
     protected:
      LocalFunctionType localFunction_;
      IndicatorType indicator_;
      BoundaryConditionsType supported_;
      std::string name_;
    };
 
 
    template<class Fct, class Indicator = EntireBoundaryIndicator,
             std::enable_if_t<(IsWrappableByConstLocalFunction<Fct>::value
                               && !GridFunction::HasRegularity<Fct, 1UL>::value
      ), int> = 0>
    constexpr auto
    neumannBoundaryModel(Fct&& values,
                         Indicator&& where = std::decay_t<Indicator>{},
                         const std::string& name = "")
    {
      return expressionClosure(NeumannBoundaryModel<Fct, Indicator>(std::forward<Fct>(values), std::forward<Indicator>(where), name));
    }
 
    template<class Fct, class Indicator = EntireBoundaryIndicator,
             std::enable_if_t<GridFunction::HasRegularity<Fct, 1UL>::value, int> = 0>
    constexpr auto
    neumannBoundaryModel(Fct&& values,
                         Indicator&& where = std::decay_t<Indicator>{},
                         const std::string& name = "")
    {
      using namespace Literals;
      return neumannBoundaryModel(gridFunction(values.gridPart(), std::forward<Fct>(values), 0_c),
                                  std::forward<Indicator>(where),
                                  name);
    }
 
    template<class Fct, class Indicator = EntireBoundaryIndicator,
             std::enable_if_t<ExpressionTraits<Fct>::isZero || ExpressionTraits<Indicator>::isZero, int> = 0>
    constexpr auto
    neumannBoundaryModel(Fct&& values,
                         Indicator&& where = std::decay_t<Indicator>{},
                         const std::string& name = "")
    {
      return zeroModel(values);
    }
 
 
 
 
  } // namespace ACFem::PDEModel::
 
  namespace ACFem {
 
    using PDEModel::neumannBoundaryModel;
 
  }
 
}  //Namespace Dune
 
 
#endif // __DUNE_ACFEM_MODELS_MODULES_NEUMANNMODEL_HH__