diracfunctional.hh

#ifndef __DUNE_ACFEM_DIRACFUNCTIONAL_HH__
#define __DUNE_ACFEM_DIRACFUNCTIONAL_HH__
 
#include "linearfunctional.hh"
#include "../../common/stringconversion.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class DiscreteFunctionSpace>
    class DiracDistribution;
 
    template<class DiscreteFunctionSpace>
    class LocalDiracDistribution;
 
    template<class DiscreteFunctionSpace>
    struct DiracDistributionTraits
      : public LinearFunctionalTraitsDefault<DiracDistribution<DiscreteFunctionSpace>,
                                             LocalDiracDistribution<DiscreteFunctionSpace> >
    {};
 
    template<class DiscreteFunctionSpace>
    class DiracDistribution
      : public DiscreteLinearFunctionalDefault<DiscreteFunctionSpace,
                                               DiracDistributionTraits<DiscreteFunctionSpace> >
    {
      typedef DiracDistribution ThisType;
      typedef
        DiscreteLinearFunctionalDefault<
          DiscreteFunctionSpace, DiracDistributionTraits<DiscreteFunctionSpace> >
        BaseType;
    public:
      typedef DiracDistributionTraits<DiscreteFunctionSpace> TraitsType;
      typedef typename TraitsType::LocalFunctionalType LocalFunctionalType;
      typedef DiscreteFunctionSpace DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
      typedef typename DiscreteFunctionSpaceType::EntityType EntityType;
      typedef typename DiscreteFunctionSpaceType::RangeFieldType FieldType;
      typedef typename DiscreteFunctionSpaceType::DomainType GlobalCoordinateType;
      typedef FieldType RangeType;
    protected:
      typedef typename GridPartType::IndexSetType IndexSetType;
      typedef typename IndexSetType::IndexType EntityIndexType;
      typedef std::vector<RangeType> CoeffsType;
      typedef typename GridPartType::template Codim<0>::GeometryType GeometryType;
      typedef typename GeometryType::LocalCoordinate LocalCoordinateType;
      typedef
        ReferenceElements<typename GeometryType::ctype, GridPartType::dimension>
        ReferenceElementsType;
 
      friend LocalDiracDistribution<DiscreteFunctionSpace>;
    public:
      DiracDistribution(const DiscreteFunctionSpaceType& space,
                        const GlobalCoordinateType& x0, unsigned component = 0)
        : BaseType(space),
          x0_(x0),
          component_(component),
          gridPart_(space_.gridPart()),
          indexSet_(gridPart_.indexSet()),
          sequence_(-1),
          found_(false),
          continuous_(space_.continuous())
      {}
 
      using BaseType::operator();
      using BaseType::coefficients;
      using BaseType::space;
      using BaseType::localFunctional;
 
      std::string name() const
      {
        std::string result;
        if (objectToString(x0_, result)) {
          return "<delta_["+result+"],...>";
        } else {
          return "<Dirac_x0,...>";
        }
      }
 
      GlobalCoordinateType& x0() { return x0_; }
      const GlobalCoordinateType& x0() const { return x0_; }
      unsigned& component() { return component_; }
      const unsigned& component() const { return component_; }
 
    protected:
      template<class Entity>
      bool doCheckEntity(const Entity& entity) const
      {
        const auto& geometry = entity.geometry();
        x0Local_ = geometry.local(x0_);
 
        if (ReferenceElementsType::general(geometry.type()).checkInside(x0Local_)) {
          return true;
        }
        return false;
      }
      template<class Entity>
      bool checkEntity(const Entity& entity) const
      {
        if (!continuous_) {
          return doCheckEntity(entity);
        }
        if (sequence_ != space_.sequence()) {
          found_ = false;
          sequence_ = space_.sequence();
        }
        if (found_) {
          return indexSet_.index(entity) == entityIndex_;
        }
#if 0
        if (doCheckEntity(entity)) {
          found_ = true;
          entityIndex_ = indexSet_.index(entity);
          return true;
        }
        return false;
#else
        // should be faster for hiearchical grids
        Fem::EntitySearch<GridPartType> entitySearch(gridPart_);
 
        const auto entityOfXZero = entitySearch(x0_);
        found_ = true; // otherwise error
        entityIndex_ = indexSet_.index(entityOfXZero);
        const auto& geometry = entityOfXZero.geometry();
        x0Local_ = geometry.local(x0_);
 
        return indexSet_.index(entity) == entityIndex_;
#endif
      }
 
    protected:
      using BaseType::space_;
      GlobalCoordinateType x0_;
      unsigned int component_;
      const GridPartType& gridPart_;
      const IndexSetType& indexSet_;
      mutable int sequence_;
      mutable EntityIndexType entityIndex_;
      mutable bool found_;
      mutable LocalCoordinateType x0Local_;
      const bool continuous_;
    };
 
    template<class DiscreteFunctionSpace>
    class LocalDiracDistribution
      : public LocalLinearFunctionalDefault<DiscreteFunctionSpace,
                                            DiracDistributionTraits<DiscreteFunctionSpace> >
    {
      typedef LocalDiracDistribution ThisType;
      typedef
      LocalLinearFunctionalDefault<DiscreteFunctionSpace,
                                   DiracDistributionTraits<DiscreteFunctionSpace> >
      BaseType;
    public:
      typedef DiracDistributionTraits<DiscreteFunctionSpace> TraitsType;
      typedef DiracDistribution<DiscreteFunctionSpace> FunctionalType;
      typedef DiscreteFunctionSpace DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::EntityType EntityType;
      typedef typename DiscreteFunctionSpaceType::RangeFieldType FieldType;
      typedef FieldType RangeType;
     public:
 
      LocalDiracDistribution(const FunctionalType& phi)
        : BaseType(phi),
          activeEntity_(false)
      {}
 
      LocalDiracDistribution(const EntityType& entity, const FunctionalType& phi)
        : BaseType(entity, phi),
          activeEntity_(false)
      {}
 
    public:
      void init(const EntityType& entity)
      {
        activeEntity_ = functional().checkEntity(entity);
        BaseType::init(entity);
      }
 
      template<class LocalFunction>
      RangeType operator()(const LocalFunction& arg) const
      {
        if (!activeEntity_) {
          return 0;
        }
        return BaseType::operator()(arg);
      }
 
      template<class LocalFunction>
      void coefficients(const RangeType& c, LocalFunction& coeffs) const
      {
        if (!activeEntity_) {
          return;
        }
        typename LocalFunction::RangeType value(0);
        value[functional().component_] = c;
        coeffs.axpy(functional().x0Local_, value);
      }
 
      template<class LocalFunction>
      void coefficients(LocalFunction& coeffs) const
      {
        BaseType::coefficients(coeffs);
      }
 
      using BaseType::operator();
      using BaseType::functional;
      using BaseType::entity;
 
    protected:
      using BaseType::entity_;
      bool activeEntity_;
    };
 
 
 
  } // namespace ACFem
 
} // namespace Dune
 
#endif // __DUNE_ACFEM_DIRACFUNCTIONAL_HH__