righthandsidefunction.hh

#ifndef DUNE_ACFEM_RIGHTHANDSIDEFUNCTION_HH
#define DUNE_ACFEM_RIGHTHANDSIDEFUNCTION_HH
 
#include "../functions/localfunctionwrapper.hh"
#include "../functions/gridfunctionexpression.hh"
#include "../models/operatorparts/operatorpartsexpression.hh"
#include "../models/modelinterface.hh"
 
namespace Dune {
  namespace ACFem {
 
    template<class OperatorParts, class Solution>
    class LocalRightHandSideAdapter
    {
     public:
      typedef OperatorPartsInterface<OperatorParts> OperatorPartsInterfaceType;
      typedef OperatorParts OperatorPartsType;
      typedef Solution SolutionType;
     protected:
      typedef typename OperatorPartsType::FunctionSpaceType OperatorPartsFunctionSpaceType;
      typedef typename SolutionType::FunctionSpaceType SolutionFunctionSpaceType;
      typedef typename SolutionType::LocalFunctionType LocalSolutionType;
 
      static_assert(std::is_same<OperatorPartsFunctionSpaceType, SolutionFunctionSpaceType>::value,
                    "OperatorParts and  \"exact\" solution must live in the same function space");
 
     public:
      typedef OperatorPartsFunctionSpaceType FunctionSpaceType; // ahem, make change in the future, maybe ...
      typedef typename SolutionType::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;
 
      LocalRightHandSideAdapter(const OperatorPartsInterfaceType& operatorParts,
                                const Fem::Function<FunctionSpaceType, SolutionType>& solution,
                                const std::string& name = "")
        : operatorParts_(operatorParts),
          solution_(solution),
          localSolution_(solution_()),
          name_(name == "" ? operatorParts_().name() + "(" + solution_().name() + ", . )" : name)
      {}
 
      LocalRightHandSideAdapter(const LocalRightHandSideAdapter& other)
        : operatorParts_(other.operatorParts_()),
          solution_(other.solution_()),
          localSolution_(solution_()),
          name_(other.name_)
      {}
 
      const std::string& name() const { return name_; }
 
      // Implement needed interface method for the adapter class
      void init(const EntityType& entity)
      {
        operatorParts_().setEntity(entity);
        localSolution_.init(entity);
      }
 
      int order () const
      {
        return localSolution_.order();
      }
 
      template<class PointType>
      void evaluate(const PointType& x, RangeType& result) const
      {
        // Now do the stuff ... we need the source term and the
        // fluxDivergence(). For the latter we need to remember the
        // sign convention we chose ages ago ... actually
        // fluxDivergence really is MINUS the divergence, i.e. we have
        // the correct sign as obtained when deriving the weak
        // formulation with integration by parts.
 
        RangeType value;
        JacobianRangeType jacobian;
        HessianRangeType hessian;
 
        localSolution_.evaluate(x, value);
        localSolution_.jacobian(x, jacobian);
        localSolution_.hessian(x, hessian);
 
        result = RangeType(0);
 
        if (OperatorPartsType::hasFlux) {
          // take the principal part ...
          operatorParts_().fluxDivergence(localSolution_.entity(), x,
                                          value, jacobian, hessian,
                                          result);
        }
 
        if (OperatorPartsType::hasSources) {
          // add the source terms
          RangeType source;
          operatorParts_().source(localSolution_.entity(), x,
                                  value, jacobian,
                                  source);
          result += source;
        }
 
      }
 
      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
        DUNE_THROW(NotImplemented, "Jacobian of an automatic RHS force density just cannot be computed.");
#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 an automatic RHS force density just cannot be computed.");
#endif
      }
 
     protected:
      ExpressionStorage<OperatorPartsType> operatorParts_;
      ExpressionStorage<SolutionType> solution_;
      mutable LocalSolutionType localSolution_;
      const std::string name_;
    };
 
 
    template<class OperatorParts, class Solution>
    LocalFunctionWrapper<LocalRightHandSideAdapter<OperatorParts, Solution>, typename Solution::GridPartType>
    rightHandSide(const OperatorPartsInterface<OperatorParts>& operatorParts,
                  const Fem::Function<typename Solution::FunctionSpaceType, Solution>& solution_,
                  const std::string& name = "")
    {
      typedef LocalRightHandSideAdapter<OperatorParts, Solution> LocalFunctionType;
      typedef LocalFunctionWrapper<LocalFunctionType, typename Solution::GridPartType> GridFunctionType;
 
      const Solution& solution(static_cast<const Solution&>(solution_));
      LocalFunctionType L_f_dot(operatorParts, solution, name);
 
      return GridFunctionType(L_f_dot.name(), L_f_dot, solution.space().gridPart(), solution.space().order());
    }
 
    template<class Model, class Solution>
    LocalFunctionWrapper<LocalRightHandSideAdapter<typename Model::OperatorPartsType, Solution>,
                         typename Solution::GridPartType>
    rightHandSide(const ModelInterface<Model>& model,
                  const Fem::Function<typename Solution::FunctionSpaceType, Solution>& solution_,
                  const std::string& name = "")
    {
      typedef LocalRightHandSideAdapter<typename Model::OperatorPartsType, Solution> LocalFunctionType;
      typedef LocalFunctionWrapper<LocalFunctionType, typename Solution::GridPartType> GridFunctionType;
 
      const Solution& solution(static_cast<const Solution&>(solution_));
      LocalFunctionType L_f_dot(model.operatorParts(), solution, name);
 
      return GridFunctionType(L_f_dot.name(), L_f_dot, solution.space().gridPart(), solution.space().order());
    }
 
 
 
  } // ACFem::
 
} // Dune::
 
 
#endif // DUNE_ACFEM_RIGHTHANDSIDEFUNCTION_HH