gradientfunction.hh

#ifndef DUNE_ACFEM_GRADIENTFUNCTION_HH
#define DUNE_ACFEM_GRADIENTFUNCTION_HH
 
#include "../functions/localfunctionwrapper.hh"
#include "../functions/gridfunctionexpression.hh"
 
namespace Dune {
  namespace ACFem {
 
    template<class GridFunction>
    class LocalGradientAdapter
    {
      typedef GridFunction GradientFunctionType;
      typedef typename GradientFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename GradientFunctionType::LocalFunctionType LocalGradientFunctionType;
     public:
      typedef typename DiscreteFunctionSpaceType::FunctionSpaceType GradientFunctionSpaceType;
      typedef typename
      Fem::ToNewDimRangeFunctionSpace<GradientFunctionSpaceType, GradientFunctionSpaceType::dimDomain>::Type
      FunctionSpaceType;
      typedef typename GridFunction::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;
 
      static_assert((int)GradientFunctionSpaceType::dimRange == 1,
                    "Gradient Adapter needs dimRange == 1 for underlying function");
 
      LocalGradientAdapter(const Fem::Function<GradientFunctionSpaceType, GradientFunctionType>& function,
                           const std::string& name = "")
        : function_(function),
          localFunction_(function_()),
          name_(name == "" ? "grad(" + function_().name() + ")" : name)
      {}
 
      LocalGradientAdapter(const LocalGradientAdapter& other)
        : function_(other.function_()),
          localFunction_(function_()),
          name_(other.name_)
      {}
 
      const std::string& name() const { return name_; }
 
      // Implement needed interface method for the adapter class
      void init(const EntityType& entity)
      {
        localFunction_.init(entity);
      }
 
      template<class PointType>
      void evaluate(const PointType& x, RangeType& ret) const
      {
        typename GradientFunctionType::JacobianRangeType myValue;
        localFunction_.jacobian(x, myValue);
 
        // Now, the gradient is always a matrix. Convert it ...
        ret = myValue[0];
      }
 
      template<class PointType>
      void jacobian(const PointType& x, JacobianRangeType& ret) const
      {
        typename GradientFunctionType::HessianRangeType myJacobian;
        localFunction_.hessian(x, myJacobian);
 
        // The Hessian is a tensor. The actual jacobian in this case
        // indeed is just the Hessian of the one and only component of
        // the original function ...
        ret = myJacobian[0];
      }
 
      // 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 gradient of a given function is NEVER implemented.");
#endif
      }
 
      int order() const
      {
        return std::max(0, localFunction_.order() - 1);
      }
 
     protected:
      ExpressionStorage<GradientFunctionType> function_;
      mutable LocalGradientFunctionType localFunction_;
      const std::string name_;
    };
 
 
 
  } // ACFem::
 
} // Dune::
 
 
#endif // DUNE_ACFEM_GRADIENTFUNCTION_HH