piecewisefunction.hh

#ifndef DUNE_ACFEM_PIECEWISEFUNCTION_HH
#define DUNE_ACFEM_PIECEWISEFUNCTION_HH
 
#include "../functions/localfunctionwrapper.hh"
#include "../functions/gridfunctionexpression.hh"
 
namespace Dune {
  namespace ACFem {
 
    template<class CharacteristicFunction, class GridFunction0, class GridFunction1>
    class LocalPiecewiseAdapter
    {
      typedef CharacteristicFunction CharacteristicFunctionType;
      typedef typename CharacteristicFunctionType::DiscreteFunctionSpaceType DiscreteCharacteristicFunctionSpaceType;
      typedef typename CharacteristicFunctionType::LocalFunctionType LocalCharacteristicFunctionType;
 
      typedef GridFunction0 PiecewiseFunctionType0;
      typedef typename PiecewiseFunctionType0::DiscreteFunctionSpaceType DiscreteFunctionSpaceType0;
      typedef typename PiecewiseFunctionType0::LocalFunctionType LocalPiecewiseFunctionType0;
      typedef GridFunction1 PiecewiseFunctionType1;
      typedef typename PiecewiseFunctionType1::DiscreteFunctionSpaceType DiscreteFunctionSpaceType1;
      typedef typename PiecewiseFunctionType1::LocalFunctionType LocalPiecewiseFunctionType1;
 
     public:
      typedef typename DiscreteFunctionSpaceType0::FunctionSpaceType FunctionSpaceType;
 
      static_assert(std::is_same<FunctionSpaceType, typename DiscreteFunctionSpaceType1::FunctionSpaceType>::value,
                    "Piecewise Adapter needs same function space for underlying functions");
      typedef typename FunctionSpaceType::ScalarFunctionSpaceType ScalarFunctionSpaceType;
      static_assert(std::is_same<ScalarFunctionSpaceType, typename DiscreteCharacteristicFunctionSpaceType::FunctionSpaceType>::value,
                    "Piecewise Adapter needs a scalar characteristic function");
 
      static_assert(std::is_base_of<Fem::HasLocalFunction, CharacteristicFunction>::value,
          "Missing HasLocalFunction property in CharacteristicFunction");
      static_assert(std::is_base_of<Fem::HasLocalFunction, GridFunction0>::value,
          "Missing HasLocalFunction property in GridFunction0");
      static_assert(std::is_base_of<Fem::HasLocalFunction, GridFunction1>::value,
          "Missing HasLocalFunction property in GridFunction1");
 
      typedef typename CharacteristicFunction::GridPartType GridPartType;
      typedef typename GridPartType::template Codim<0>::EntityType EntityType;
 
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
      LocalPiecewiseAdapter(const Fem::Function<ScalarFunctionSpaceType,
                            CharacteristicFunctionType>& characteristic,
                            const Fem::Function<FunctionSpaceType,
                            PiecewiseFunctionType0>& function0,
                            const Fem::Function<FunctionSpaceType,
                            PiecewiseFunctionType1>& function1,
                            const std::string& name = "")
        : characteristicFunction_(characteristic),
          localCharacteristicFunction_(characteristicFunction_()),
          function0_(function0),
          localFunction0_(function0_()),
          function1_(function1),
          localFunction1_(function1_()),
          name_(name == "" ? "(" + characteristicFunction_().name() + "*" + function1_().name() + "+" + "(1-" + characteristicFunction_().name() + ")*" + function0_().name() +")" : name)
      {}
 
      LocalPiecewiseAdapter(const LocalPiecewiseAdapter& other)
        : characteristicFunction_(other.characteristicFunction_()),
          localCharacteristicFunction_(characteristicFunction_()),
          function0_(other.function0_()),
          localFunction0_(function0_()),
          function1_(other.function1_()),
          localFunction1_(function1_()),
          name_(other.name_)
      {}
 
      const std::string& name() const { return name_; }
 
      // Implement needed interface method for the adapter class
      void init(const EntityType& entity)
      {
        localCharacteristicFunction_.init(entity);
        localFunction0_.init(entity);
        localFunction1_.init(entity);
      }
 
      template<class PointType>
      void evaluate(const PointType& x, RangeType& ret) const
      {
        typename ScalarFunctionSpaceType::RangeType indicator;
        localCharacteristicFunction_.evaluate(x, indicator);
        if (indicator > 0.5)
          localFunction1_.evaluate(x, ret);
        else
          localFunction0_.evaluate(x, ret);
      }
 
      template<class PointType>
      void jacobian(const PointType& x, JacobianRangeType& ret) const
      {
#if __DUNE_ACFEM_MAKE_CHECK__
        // in order not to disturb the test-suite, which greedily takes hessians from anything
        ret = 0;
#else
        typename ScalarFunctionSpaceType::RangeType indicator;
        localCharacteristicFunction_.evaluate(x, indicator);
        if (indicator > 0.5)
          localFunction1_.jacobian(x, ret);
        else
          localFunction0_.jacobian(x, ret);//DUNE_THROW(NotImplemented, "Jacobian of the piecewise of a given function is a singularity.");
#endif
      }
 
      // hessian of local function
      template<class PointType>
      void hessian(const PointType& x, HessianRangeType& ret) const
      {
#if __DUNE_ACFEM_MAKE_CHECK__
        // in order not to disturb the test-suite, which greedily takes hessians from anything
        ret = 0;
#else
        DUNE_THROW(NotImplemented, "Hessian of the piecewise of a given function is NEVER implemented.");
#endif
      }
 
      int order() const
      {
        // TODO: make use of characteristic function here.
        //       what about boundary elements?
        return std::max(localFunction0_.order(), localFunction1_.order());
      }
 
     protected:
      ExpressionStorage<CharacteristicFunctionType> characteristicFunction_;
      mutable LocalCharacteristicFunctionType localCharacteristicFunction_;
      ExpressionStorage<PiecewiseFunctionType0> function0_;
      mutable LocalPiecewiseFunctionType0 localFunction0_;
      ExpressionStorage<PiecewiseFunctionType1> function1_;
      mutable LocalPiecewiseFunctionType1 localFunction1_;
      const std::string name_;
    };
 
 
 
  } // ACFem::
 
} // Dune::
 
 
#endif // DUNE_ACFEM_PIECEWISEFUNCTION_HH