divergencefunction.hh

#ifndef DUNE_ACFEM_DIVERGENCEFUNCTION_HH
#define DUNE_ACFEM_DIVERGENCEFUNCTION_HH
 
#include "../functions/localfunctionwrapper.hh"
#include "../functions/gridfunctionexpression.hh"
 
namespace Dune {
  namespace ACFem {
 
    template<class GridFunction>
    class LocalDivergenceAdapter
    {
      typedef GridFunction DivergenceFunctionType;
      typedef typename DivergenceFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename DivergenceFunctionType::LocalFunctionType LocalDivergenceFunctionType;
     public:
      typedef typename DiscreteFunctionSpaceType::FunctionSpaceType DivergenceFunctionSpaceType;
      typedef typename DivergenceFunctionSpaceType::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;
 
      static_assert((int)DivergenceFunctionSpaceType::dimRange == (int)FunctionSpaceType::dimDomain,
                    "Divergence Adapter needs dimRange == dimDomain for underlying function");
 
      LocalDivergenceAdapter(const Fem::Function<DivergenceFunctionSpaceType, DivergenceFunctionType>& function,
                             const std::string& name = "")
        : function_(function),
          localFunction_(function_()),
          name_(name == "" ? "div(" + function_().name() + ")" : name)
      {}
 
      LocalDivergenceAdapter(const LocalDivergenceAdapter& 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 DivergenceFunctionType::JacobianRangeType divJacobian;
        localFunction_.jacobian(x, divJacobian);
 
        ret = 0;
        for (int i = 0; i < FunctionSpaceType::dimDomain; ++i) {
          ret[0] += divJacobian[i][i];
        }
      }
 
      template<class PointType>
      void jacobian(const PointType& x, JacobianRangeType& ret) const
      {
        typename DivergenceFunctionType::HessianRangeType divHessian;
        localFunction_.hessian(x, divHessian);
 
        ret = 0;
        for (int i = 0; i < FunctionSpaceType::dimDomain; ++i) {
          for (int j = 0; j < FunctionSpaceType::dimDomain; ++j) {
            ret[0][i] += divHessian[0][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 divergence of a given function is NEVER implemented.");
#endif
      }
 
     protected:
      ExpressionStorage<DivergenceFunctionType> function_;
      mutable LocalDivergenceFunctionType localFunction_;
      const std::string name_;
    };
 
 
 
  } // ACFem::
 
} // Dune::
 
 
#endif // DUNE_ACFEM_DIVERGENCEFUNCTION_HH