probleminterface.hh

#ifndef _DUNE_ACFEM_POISSON_PROBLEMINTERFACE_HH_
#define _DUNE_ACFEM_POISSON_PROBLEMINTERFACE_HH_
 
#include <cassert>
#include <cmath>
 
#include <dune/common/exceptions.hh>
#include <dune/fem/function/common/function.hh>
#include "timeview.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template <class FunctionSpace>
    class ProblemInterface
    {
     public:
      // type of function space
      typedef FunctionSpace  FunctionSpaceType;
 
      enum { dimRange  = FunctionSpaceType::dimRange  };
      enum { dimDomain = FunctionSpaceType::dimDomain };
 
      typedef typename FunctionSpaceType::RangeFieldType   RangeFieldType;
 
      typedef typename FunctionSpaceType::RangeType   RangeType;
      typedef typename FunctionSpaceType::DomainType  DomainType;
 
      typedef typename FunctionSpaceType::JacobianRangeType  JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
      enum OperatorParts {
        SecondOrderCoefficient,
        FirstOrderCoefficient,
        ZeroOrderCoefficient,
        RightHandSide,
        DirichletBoundary,
        NeumannBoundary,
        RobinBoundary,
        ExactSolution,
      };
      typedef enum OperatorParts OperatorPartsType;
 
      virtual ~ProblemInterface() {}
 
      virtual bool has(OperatorPartsType what) const
      {
        switch (what) {
        case SecondOrderCoefficient: return true;
        case FirstOrderCoefficient:  return false;
        case ZeroOrderCoefficient:   return false;
        case RightHandSide:          return true;
        case DirichletBoundary:      return true;
        case NeumannBoundary:        return false;
        case RobinBoundary:          return false;
        case ExactSolution:          return true;
        }
        return false;
      }
 
      virtual void f(const DomainType& x,
                     RangeType& value) const
      {
        value = 0;
      }
 
      virtual void u(const DomainType& x,
                     RangeType& value) const
      {
        value = 0;
      }
 
      virtual void uJacobian(const DomainType& x,
                             JacobianRangeType& value) const
      {
        value = 0;
      }
 
      virtual bool isDirichletSegment(const int bndId, const DomainType& center) const {
        return true;
      }
 
      virtual bool isNeumannSegment(const int bndId, const DomainType& center) const {
        return false;
      }
 
      virtual bool isRobinSegment(const int bndId, const DomainType& center) const {
        return false;
      }
 
      virtual void dirichletData(const DomainType& x,
                                 RangeType& value) const
      {
        u(x, value);
      }
 
      virtual void neumannData(const DomainType& x,
                               RangeType& value) const
      {
        value = 0;
      }
 
      virtual void robinData(const DomainType& x,
                             const RangeType& u,
                             RangeType& value) const
      {
        value = 0;
      }
 
 
 
      virtual void secondOrderCoefficient(const DomainType& x,
                                          const JacobianRangeType& gradient,
                                          JacobianRangeType& result) const
      {
        // default implementation is Laplace
        result = gradient ;
      }
 
      virtual void secondOrderCoefficient(const DomainType &x,
                                          const JacobianRangeType &Du,
                                          const HessianRangeType &D2u,
                                          RangeType &result) const
      {
        // default implementation is Laplace
        for (int i = 0; i < dimRange; ++i) {
          result[i] = 0;
          for (int j = 0; j < dimDomain; ++j) {
            result[i] -= D2u[i][j][j];
          }
        }
      }
 
      virtual void firstOrderCoefficient(const DomainType& x,
                                         const RangeType& u,
                                         const JacobianRangeType& Du,
                                         RangeType& result) const
      {
        result = 0;
      }
 
 
      virtual void zeroOrderCoefficient(const DomainType& x,
                                        const RangeType& u,
                                        RangeType& result) const
      {
        result = 0;
      }
 
      enum FunctionId { exact, rhs, dirichlet, neumann };
 
      template <FunctionId id>
      class FunctionWrapper : public Dune::Fem::Function<FunctionSpaceType, FunctionWrapper<id> >
      {
        const ProblemInterface& impl_;
 
       public:
        FunctionWrapper(const ProblemInterface& impl)
          : impl_(impl)
        {
          //std::cout << "ctor impl: " << impl_ << " " << &impl_ << " " << this << std::endl;
        }
 
        void evaluate(const DomainType& x, RangeType& ret) const
        {
          //std::cout << "eval impl: " << impl_ << " " << &impl_ << " " << this << std::endl;
          switch (id) {
          case exact: // call exact solution of implementation
            impl_.u(x, ret);
            break;
          case rhs: // call right hand side of implementation
            impl_.f(x, ret);
            break;
          case dirichlet: // call dirichlet boudary data of implementation
            impl_.dirichletData(x, ret);
            break;
          case neumann: // call neumann boundary data of implementation
            impl_.neumannData(x, ret);
            break;
          default:
            DUNE_THROW(Dune::NotImplemented,"FunctionId not implemented");
            break;
          }
        }
 
        void jacobian(const DomainType& x, JacobianRangeType& jac) const
        {
          switch (id) {
          case exact:
            impl_.uJacobian(x, jac);
            break;
          default:
            DUNE_THROW(Dune::NotImplemented,"Jacobian for FunctionId not implemented");
            break;
          }
        }
      };
 
      typedef FunctionWrapper<exact> ExactSolutionType;
      // return Fem::Function for exact solution
      ExactSolutionType exactSolution() const
      {
        return ExactSolutionType(*this);
      }
 
      typedef FunctionWrapper<rhs> RightHandSideType;
      // return Fem::Function for right hand side
      RightHandSideType rightHandSide() const
      {
        return RightHandSideType(*this);
      }
 
      typedef FunctionWrapper<dirichlet> DirichletBoundaryType;
      // return Fem::Function for Dirichlet boundary values
      DirichletBoundaryType dirichletBoundary() const
      {
        return DirichletBoundaryType(*this);
      }
 
      typedef FunctionWrapper<neumann> NeumannBoundaryType;
      // return Fem::Function for Neumann boundary values
      NeumannBoundaryType neumannBoundary() const
      {
        return NeumannBoundaryType(*this);
      }
    };
 
    template <class FunctionSpace, class TimeProvider>
    class TransientProblemInterface
      : public ProblemInterface<FunctionSpace>
    {
      typedef ProblemInterface<FunctionSpace> BaseType;
      typedef TransientProblemInterface<FunctionSpace, TimeProvider> ThisType;
     public:
      typedef FunctionSpace FunctionSpaceType;
      typedef TimeProvider TimeProviderType;
      typedef Dune::ACFem::TimeView<TimeProviderType> TimeViewType;
 
      TransientProblemInterface(const TimeProviderType& timeProvider, double theta = 0.0)
        : timeView_(timeProvider, theta)
      {}
 
      typedef typename BaseType::ExactSolutionType InitialValueType;
 
      // return Fem::Function for initial
      InitialValueType initialValue() const
      {
        return InitialValueType(*this);
      }
 
      double time() const
      {
        return timeView_.time();
      }
 
      double deltaT() const
      {
        return timeView_.deltaT() ;
      }
 
      const TimeViewType& timeView() const
      {
        return timeView_;
      }
 
      TimeViewType& timeView()
      {
        return timeView_;
      }
     protected:
      TimeViewType timeView_;
    };
 
 
 
  } // ACFem::
 
} // Dune::
 
#endif // #ifndef ELLIPTC_PROBLEMINTERFACE_HH