maxfunction.hh

#ifndef DUNE_ACFEM_MAXFUNCTION_HH
#define DUNE_ACFEM_MAXFUNCTION_HH
 
#include "constantfunction.hh"
 
namespace Dune {
  namespace ACFem {
 
    template<class GridFunction1, class GridFunction2>
    class LocalMaxAdapter
    {
      static_assert(std::is_same<typename GridFunction1::FunctionSpaceType,
          typename GridFunction2::FunctionSpaceType>::value, "Functions have different FunctionSpaceType");
 
      typedef typename GridFunction1::LocalFunctionType LocalFunctionType1;
      typedef typename GridFunction2::LocalFunctionType LocalFunctionType2;
 
     public:
      typedef typename GridFunction1::FunctionSpaceType FunctionSpaceType;
      typedef typename GridFunction1::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;
      enum{dimRange = FunctionSpaceType::dimRange};
 
      LocalMaxAdapter(const Fem::Function<FunctionSpaceType, GridFunction1>& f1,
                      const Fem::Function<FunctionSpaceType, GridFunction2>& f2,
                      const std::string& name = "")
        : f1_(f1), f2_(f2),
          locf1_(f1_()), locf2_(f2_()),
          name_(name == "" ? "max{" + f1_().name() + "," + f2_().name() + "}" : name)
      {}
 
      LocalMaxAdapter(const LocalMaxAdapter& other)
        : f1_(other.f1_()), f2_(other.f2_()),
          locf1_(f1_()), locf2_(f2_()),
          name_(other.name_)
      {}
 
      const std::string& name() const
      {
        return name_;
      }
 
      void init(const EntityType& entity)
      {
        locf1_.init(entity);
        locf2_.init(entity);
      }
 
      template<class PointType>
      void evaluate(const PointType& x, RangeType& ret) const
      {
        RangeType value2;
 
        locf1_.evaluate(x, ret);
        locf2_.evaluate(x, value2);
        for(int i = 0; i < dimRange; ++ i)
          if (ret[i] < value2[i])
            ret[i] = value2[i];
      }
 
      template<class PointType>
      void jacobian(const PointType& x, JacobianRangeType& ret) const
      {
        RangeType value1;
        RangeType value2;
        locf1_.evaluate(x, value1);
        locf2_.evaluate(x, value2);
        JacobianRangeType jacValue2;
        locf1_.jacobian(x, ret);
        locf2_.jacobian(x, jacValue2);
 
        for(int i = 0; i < dimRange; ++i)
          if (value1[i] < value2[i])
            ret[i] = jacValue2[i];
      }
 
      template<class PointType>
      void hessian(const PointType& x, HessianRangeType& ret) const
      {
        RangeType value1;
        RangeType value2;
        locf1_.evaluate(x, value1);
        locf2_.evaluate(x, value2);
        HessianRangeType hesValue2;
        locf1_.hessian(x, ret);
        locf2_.hessian(x, hesValue2);
 
        for(int i = 0; i < dimRange; ++i)
          if (value1[i] < value2[i])
            ret[i] = hesValue2[i];
      }
 
      int order() const
      {
        return std::max(locf1_.order(), locf2_.order());
      }
 
     protected:
      ExpressionStorage<GridFunction1> f1_;
      ExpressionStorage<GridFunction2> f2_;
      mutable LocalFunctionType1 locf1_;
      mutable LocalFunctionType2 locf2_;
      const std::string name_;
    };
 
 
    template<class GridFunction1, class GridFunction2>
    LocalFunctionWrapper<LocalMaxAdapter<GridFunction1, GridFunction2>, typename GridFunction1::GridPartType>
    max(const Fem::Function<typename GridFunction1::FunctionSpaceType, GridFunction1>& f1,
        const Fem::Function<typename GridFunction2::FunctionSpaceType, GridFunction2>& f2,
         const std::string& name = "")
    {
      const auto& f1_ = static_cast<const GridFunction1&>(f1);
      const auto& f2_ = static_cast<const GridFunction2&>(f2);
      LocalMaxAdapter<GridFunction1, GridFunction2> local(f1_, f2_, name);
 
      return LocalFunctionWrapper<LocalMaxAdapter<GridFunction1, GridFunction2>, typename GridFunction1::GridPartType>
        (local.name(), local, f1_.space().gridPart(), std::max(f1_.space().order(), f2_.space().order()));
    }
 
    template<class GridFunction1>
    LocalFunctionWrapper<LocalMaxAdapter<GridFunction1, ConstantGridFunction< typename GridFunction1::FunctionSpaceType, typename GridFunction1::GridPartType> >, typename GridFunction1::GridPartType>
    max(const Fem::Function<typename GridFunction1::FunctionSpaceType, GridFunction1>& f1,
        const typename GridFunction1::FunctionSpaceType::RangeType value2,
         const std::string& name = "")
    {
      typedef ConstantGridFunction< typename GridFunction1::FunctionSpaceType, typename GridFunction1::GridPartType> GridFunction2;
      const auto& f1_ = static_cast<const GridFunction1&>(f1);
      const GridFunction2 f2_(value2 , f1_.space().gridPart() );
      LocalMaxAdapter<GridFunction1, GridFunction2> local(f1_, f2_, name);
 
      return LocalFunctionWrapper<LocalMaxAdapter<GridFunction1, GridFunction2>, typename GridFunction1::GridPartType>
        (local.name(), local, f1_.space().gridPart(), std::max(f1_.space().order(), f2_.space().order()));
    }
 
 
 
  } // ACFem::
 
} // Dune::
 
 
#endif // DUNE_ACFEM_MAXFUNCTION_HH