diracfunctional.hh

#ifndef __DUNE_ACFEM_OPERATORS_FUNCTIONALS_MODULES_DIRACFUNCTIONAL_HH__
#define __DUNE_ACFEM_OPERATORS_FUNCTIONALS_MODULES_DIRACFUNCTIONAL_HH__
 
#include <dune/fem/gridpart/common/entitysearch.hh>
 
#include "../../../common/tostring.hh"
#include "../../../common/literals.hh"
#include "../../../expressions/terminal.hh"
#include "../linearfunctional.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    namespace LinearFunctional {
 
      using namespace Literals;
 
      template<class DiscreteFunctionSpace>
      class DiracDistribution
        : public LinearFunctionalBase<DiscreteFunctionSpace>
        , Expressions::SelfExpression<DiracDistribution<DiscreteFunctionSpace> >
      {
        using BaseType = LinearFunctionalBase<DiscreteFunctionSpace>;
       public:
        using typename BaseType::DiscreteFunctionSpaceType;
        using typename BaseType::RangeFieldType;
        using BaseType::space;
        using RangeType = typename DiscreteFunctionSpaceType::RangeType;
        using DomainType = typename DiscreteFunctionSpaceType::DomainType;
        using GridPartType = typename DiscreteFunctionSpaceType::GridPartType;
 
        template<class X0Arg, std::enable_if_t<std::is_constructible<DomainType, X0Arg>::value, int> = 0>
        DiracDistribution(X0Arg&& x0, const DiscreteFunctionSpaceType& space)
          : BaseType(space)
          , x0_(std::forward<X0Arg>(x0))
        {}
 
        template<class T>
        DiracDistribution(const std::initializer_list<T> &values, const DiscreteFunctionSpaceType& space)
          : BaseType(space)
          , x0_(values)
        {}
 
        template<class FieldArg, class DiscreteFunction,
                 std::enable_if_t<std::is_constructible<RangeFieldType, FieldArg>::value, int> = 0>
        void coefficients(FieldArg&& s, DiscreteFunction& values) const
        {
          if constexpr (!ExpressionTraits<FieldArg>::isZero) {
            try {
              const auto entity = Fem::EntitySearch<GridPartType>(space().gridPart())(x0_);
              values.localFunction(entity).axpy(entity.geometry().local(x0_), RangeType(s));
            } catch (const Dune::Exception& e) {
              // That's ok. The point is just not located inside the mesh.
            }
          }
        }
 
        template<class DiscreteFunction>
        void coefficients(DiscreteFunction& df) const
        {
          coefficients(1_f, df);
        }
 
        template<class DiscreteFunction>
        auto operator()(const DiscreteFunction& df) const
        {
          return inner(*this, df);
        }
 
        std::string name() const
        {
          return "<delta_["+toString(x0_)+"],...>";
        }
 
        const auto& x0() const& { return x0_; }
        auto& x0() & { return x0_; }
        DomainType x0() && { return x0_; }
 
       protected:
        DomainType x0_;
      };
 
      template<class X0Arg, class DiscreteFunctionSpace,
               std::enable_if_t<std::is_constructible<typename DiscreteFunctionSpace::DomainType, X0Arg>::value, int> = 0>
      auto diracFunctional(X0Arg&& x0, const DiscreteFunctionSpace& space)
      {
        return expressionClosure(DiracDistribution<DiscreteFunctionSpace>(std::forward<X0Arg>(x0), space));
      }
 
      template<class T, class DiscreteFunctionSpace>
      auto diracFunctional(const std::initializer_list<T>& x0, const DiscreteFunctionSpace& space)
      {
        return DiracDistribution<DiscreteFunctionSpace>(x0, space);
      }
 
 
 
    } // LinearFunctional::
 
    using LinearFunctional::diracFunctional;
 
  } // namespace ACFem
 
} // namespace Dune
 
#endif // __DUNE_ACFEM_OPERATORS_FUNCTIONALS_MODULES_DIRACFUNCTIONAL_HH__