spaceoperatorif.hh

#ifndef DUNE_FEM_SPACEOPERATORIF_HH
#define DUNE_FEM_SPACEOPERATORIF_HH
 
//- system includes
#include <cassert>
#include <cstdlib>
#include <limits>
#include <utility>
 
//-Dune fem includes
#include <dune/fem/operator/common/automaticdifferenceoperator.hh>
#include <dune/fem/operator/common/objpointer.hh>
#include <dune/fem/function/adaptivefunction.hh>
 
namespace Dune
{
 
  namespace Fem
  {
    template< class DiscreteFunction,
              class JacobianOperator = Fem::AutomaticDifferenceOperator< DiscreteFunction >  >
    class SpaceOperatorInterface
    : public JacobianOperator
    {
      typedef SpaceOperatorInterface< DiscreteFunction, JacobianOperator > ThisType;
      typedef Fem::Operator< DiscreteFunction > BaseType;
 
    public:
      typedef DiscreteFunction DestinationType;
 
      typedef typename DestinationType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef DiscreteFunctionSpaceType SpaceType;
 
      virtual ~SpaceOperatorInterface() {}
 
      using BaseType::operator ();
 
      virtual const DiscreteFunctionSpaceType& space() const = 0;
 
      virtual int size () const { return space().size(); }
 
      virtual void initializeTimeStepSize ( const DestinationType& U0 ) const;
 
      virtual bool hasLimiter () const { return false ; }
 
      virtual void setTime ( const double time ) {}
 
      virtual double timeStepEstimate () const
      {
        return std::numeric_limits< double >::max();
      }
 
      virtual void limit (const DestinationType& arg, DestinationType& dest) const
      {
        // if this method is not overloaded then hasLimiter should return false
        assert( ! hasLimiter () );
        // default operation is copy arg to dest
        dest.assign( arg );
      }
 
      virtual void applyLimiter ( DestinationType& U ) const
      {
        if( hasLimiter() )
        {
          if( ! uTmp_ )
            uTmp_.reset( new DestinationType( "SpaceOpIF::uTmp_", space() ) );
 
          assert( uTmp_ );
          uTmp_->assign( U );
          limit( *uTmp_, U );
        }
      }
 
      virtual const DestinationType* destination() const { return nullptr; }
 
    protected:
      mutable std::unique_ptr< DestinationType > uTmp_;
    };
 
    template <class OperatorType>
    class SpaceOperatorStorage
    : public ObjPointerStorage
    {
      SpaceOperatorStorage(const SpaceOperatorStorage& org);
      SpaceOperatorStorage& operator = (const SpaceOperatorStorage& org);
 
    protected:
      // operator storage
      mutable OperatorType* op_;
      // model storage
      ObjPointerStorage* model_;
 
    public:
      SpaceOperatorStorage(OperatorType * op)
        : op_(op), model_(0)
      {}
 
      SpaceOperatorStorage(OperatorType * op, ObjPointerStorage* model)
        : op_(op), model_(model)
      {}
 
      ~SpaceOperatorStorage()
      {
        // delete operator before destructor of base class is called
        delete op_; op_ = 0;
        delete model_; model_ = 0;
      }
 
      OperatorType& pass() const
      {
        assert( op_ );
        return (*op_);
      }
    };
 
    template <class OperatorType>
    class SpaceOperatorPtr
    : public SpaceOperatorStorage< OperatorType >,
      public SpaceOperatorInterface<typename OperatorType::DestinationType>
    {
      typedef SpaceOperatorStorage< OperatorType > BaseType;
 
                        protected:
      // use pass method of base
      using BaseType :: pass;
                        private:
 
      typedef typename OperatorType::DestinationType DestinationType;
 
      typedef typename DestinationType :: DiscreteFunctionSpaceType SpaceType;
 
      SpaceOperatorPtr(const SpaceOperatorPtr& org);
      SpaceOperatorPtr& operator = (const SpaceOperatorPtr& org);
 
    public:
      SpaceOperatorPtr(OperatorType * op)
        : BaseType(op)
      {}
 
      SpaceOperatorPtr(OperatorType * op, ObjPointerStorage* model)
        : BaseType(op,model)
      {}
 
      virtual ~SpaceOperatorPtr() {}
 
      virtual void operator () (const DestinationType& arg, DestinationType& dest) const
      {
        // this method should not be called
        assert(false);
        abort();
      }
 
      const SpaceType& space() const { return pass().space(); }
 
      void setTime(const double time) { pass().setTime(time); }
 
      double timeStepEstimate () const { return pass().timeStepEstimate(); }
 
      const DestinationType* destination() const
      {
        pass().allocateLocalMemory();
        return & (pass().destination());
      }
    };
 
    template <class OperatorType>
    class SpaceOperatorWrapper
    : public SpaceOperatorPtr< OperatorType >
    {
      typedef SpaceOperatorPtr< OperatorType > BaseType;
 
      // use pass method of base
      using BaseType :: pass;
 
      SpaceOperatorWrapper(const SpaceOperatorWrapper& org);
      SpaceOperatorWrapper& operator = (const SpaceOperatorWrapper& org);
    public:
      typedef typename OperatorType::DestinationType DestinationType;
      typedef typename DestinationType :: DiscreteFunctionSpaceType SpaceType;
 
      SpaceOperatorWrapper(OperatorType * op)
        : BaseType(op)
      {}
 
      SpaceOperatorWrapper(OperatorType * op, ObjPointerStorage* model)
        : BaseType(op,model)
      {}
 
      void operator () (const DestinationType& arg, DestinationType& dest) const
      {
        pass()(arg,dest);
      }
    };
 
 
    // Implementation of SpaceOperatorInterface
    // ----------------------------------------
 
    template< class DiscreteFunction, class JacobianOperator >
    inline void SpaceOperatorInterface< DiscreteFunction, JacobianOperator >
      ::initializeTimeStepSize ( const DestinationType &U0 ) const
    {
      // create temporary variable
      DestinationType tmp( U0 );
      // call operator
      (*this)( U0, tmp );
    }
 
  } // namespace Fem
 
} // namespace Dune
 
#endif // #ifndef DUNE_FEM_SPACEOPERATORIF_HH