normfunction.hh

#ifndef DUNE_ACFEM_NORMFUNCTION_HH
#define DUNE_ACFEM_NORMFUNCTION_HH
 
 
namespace Dune {
  namespace ACFem {
 
    // TODO: implement
    template<class GridFunction, unsigned int p>
    class LocalNormAdapter;
 
 
    template<class GridFunction>
    class LocalNormAdapter<GridFunction, 2>
    {
      typedef GridFunction NormFunctionType;
      typedef typename NormFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename NormFunctionType::LocalFunctionType LocalNormFunctionType;
     public:
      typedef typename DiscreteFunctionSpaceType::FunctionSpaceType NormFunctionSpaceType;
      typedef typename NormFunctionSpaceType::ScalarFunctionSpaceType 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;
 
      LocalNormAdapter(const Fem::Function<NormFunctionSpaceType, NormFunctionType>& function,
                       const std::string& name = "")
        : function_(function),
          localFunction_(function_()),
          name_(name == "" ? "|" + function_().name() + "|_2" : name)
      {}
 
      LocalNormAdapter(const LocalNormAdapter& other)
        : function_(other.function_()),
          localFunction_(function_()),
          name_(other.name_)
      {}
 
      const std::string& name() const { return name_; }
 
      void init(const EntityType& entity)
      {
        localFunction_.init(entity);
      }
 
      template<class PointType>
      void evaluate(const PointType& x, RangeType& ret) const
      {
        typename NormFunctionType::RangeType value;
        localFunction_.evaluate(x, value);
        ret = value.two_norm();
      }
 
      template<class PointType>
      void jacobian(const PointType& x, JacobianRangeType& ret) const
      {
        typename NormFunctionType::RangeType value;
        typename NormFunctionType::JacobianRangeType jacobian;
        localFunction_.evaluate(x, value);
        localFunction_.jacobian(x, jacobian);
        const auto norm = value.two_norm();
 
        ret = 0;
        for (int j = 0; j < NormFunctionSpaceType::dimDomain; ++j) {
          for (int i = 0; i < NormFunctionSpaceType::dimRange; ++i) {
            ret[0][j] += (value[i] / norm) * jacobian[i][j];
          }
        }
      }
 
      // 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 norm of a given function is not implemented.");
#endif
      }
 
      int order() const
      {
        return localFunction_.order() == 0 ? 0 : 111;
      }
 
     protected:
      ExpressionStorage<NormFunctionType> function_;
      mutable LocalNormFunctionType localFunction_;
      const std::string name_;
    };
 
 
    template<class GridFunction>
    LocalFunctionWrapper<LocalNormAdapter<GridFunction, 2>, typename GridFunction::GridPartType>
    twoNorm(const Fem::Function<typename GridFunction::FunctionSpaceType, GridFunction>& f,
         const std::string& name = "")
    {
      const auto& f_ = static_cast<const GridFunction&>(f);
      LocalNormAdapter<GridFunction, 2> local(f_, name);
 
      return LocalFunctionWrapper<LocalNormAdapter<GridFunction, 2>, typename GridFunction::GridPartType>
        (local.name(), local, f_.space().gridPart(), f_.space().order() == 0 ? 0 : 111);
    }
 
 
  } // ACFem::
 
} // Dune::
 
 
#endif // DUNE_ACFEM_NORMFUNCTION_HH