spacetimehandler_base.h

Go to the documentation of this file.

#ifndef SPACE_TIME_HANDLER_BASE_H_
#define SPACE_TIME_HANDLER_BASE_H_

#include <deal.II/lac/vector.h>
#include <deal.II/lac/block_vector_base.h>
#include <deal.II/lac/block_vector.h>

#include <vector>
#include <iostream>
#include <sstream>
#include <limits>

#include "timedofhandler.h"
#include "timeiterator.h"
#include "refinementcontainer.h"
#include "dopetypes.h"

namespace DOpE
{
template<typename VECTOR>//has to be a class template because you can not have virtual member function templates.
class SpaceTimeHandlerBase
{
  public:

    SpaceTimeHandlerBase(DOpEtypes::ControlType control_type = DOpEtypes::stationary) : control_type_(control_type)
    {
      time_triangulation_ = NULL;
      state_ticket_ = 1;
      control_ticket_ = 1;
    }

    SpaceTimeHandlerBase(dealii::Triangulation<1> & times, DOpEtypes::ControlType type = DOpEtypes::stationary) :
      tdfh_(times), interval_(tdfh_.first_interval()), control_type_(type)
      {
        time_triangulation_ = &times;
        state_ticket_ = 1;
        control_ticket_ = 1;
      }

    SpaceTimeHandlerBase(dealii::Triangulation<1> & times,
        const dealii::FiniteElement<1>& fe,
        DOpEtypes::ControlType type = DOpEtypes::stationary) :
          tdfh_(times, fe), interval_(tdfh_.first_interval()), control_type_(type)
    {
      time_triangulation_ = &times;
      state_ticket_ = 1;
      control_ticket_ = 1;
    }


    virtual ~SpaceTimeHandlerBase()
    {
      tdfh_.clear();
    }

    void ReInitTime()
    {
      IncrementStateTicket();
      tdfh_.distribute_dofs();
      //FIXME When we have temporal discretization also for control and constraint,
      //one has to increment here the control_ticket_!
    }

    unsigned int GetMaxTimePoint() const
    {
      return (tdfh_.GetNbrOfDoFs()-1);//because we start counting at 0.
    }

    unsigned int GetNbrOfIntervals() const
    {
      return tdfh_.GetNbrOfIntervals();
    }


    void SetInterval(const TimeIterator& it)
    {
      interval_ = it;
    }


    const TimeIterator& GetInterval() const
    {
      return interval_;
    }


    double GetTime(unsigned int time_point) const
    {

      return tdfh_.GetTime(time_point);
    }


    const TimeDoFHandler& GetTimeDoFHandler() const
    {
      return tdfh_;
    }


    const std::vector<double>& GetTimes() const
    {
      return tdfh_.GetTimes();
    }

    void
    GetTimes(const TimeIterator& interval, std::vector<double>& local_times) const
    {
      return tdfh_.GetTimes(interval, local_times);
    }


    bool IsValidStateTicket(unsigned int& ticket) const
    {
      bool ret = (ticket == state_ticket_);
      ticket = state_ticket_;
      return ret;
    }

    bool IsValidControlTicket(unsigned int& ticket) const
    {
      bool ret = (ticket == control_ticket_);
      ticket = control_ticket_;
      return ret;
    }


    DOpEtypes::ControlType GetControlType() const
    {
      return control_type_;
    }


    virtual void
      InterpolateControl(VECTOR& /*result*/, const std::vector<VECTOR*> &/*local_vectors*/,
          double /*t*/, const TimeIterator&/*interval*/) const { abort(); }
    virtual void
      InterpolateState(VECTOR& /*result*/, const std::vector<VECTOR*> &/*local_vectors*/,
                       double /*t*/, const TimeIterator&/*interval*/) const = 0;
    virtual void InterpolateConstraint(VECTOR& /*result*/, const std::vector<VECTOR*> &/*local_vectors*/,
        double /*t*/, const TimeIterator&/*interval*/) const { abort(); }

    virtual unsigned int GetControlNDoFs(int /*time_point*/ = -1) const { abort(); }
    virtual unsigned int GetStateNDoFs(int time_point = -1) const = 0;
    virtual unsigned int GetConstraintNDoFs(std::string /*name*/) const { abort(); return 0; }
    virtual const std::vector<unsigned int>& GetControlDoFsPerBlock(int /*time_point*/ = -1) const  { abort(); }
    virtual const std::vector<unsigned int>& GetStateDoFsPerBlock(int time_point = -1) const = 0;
    virtual const std::vector<unsigned int>& GetConstraintDoFsPerBlock(std::string /*name*/) const { abort(); }
    virtual unsigned int  GetNGlobalConstraints() const { abort(); }
    virtual unsigned int  GetNLocalConstraints() const { abort(); }
    
    double GetStepSize() const
    {
      return interval_.get_k();
    }
    double GetNextStepSize() const
    {
      assert(interval_!= tdfh_.last_interval());
      double k = (++interval_).get_k();
      --interval_;
      return k;
    }
    double GetPreviousStepSize() const
    {
      assert(interval_!= tdfh_.first_interval());
      double k = (--interval_).get_k();
      ++interval_;
      return k;
    }


    virtual unsigned int NewTimePointToOldTimePoint(unsigned int t) const = 0;


    virtual void SpatialMeshTransferControl(const VECTOR& /*old_values*/, VECTOR& /*new_values*/) const { abort(); }
    
    virtual void SpatialMeshTransferState(const VECTOR& /*old_values*/, VECTOR& /*new_values*/) const { abort(); }

    /******************************************************/
        void
          RefineTime(DOpEtypes::RefinementType /*ref_type*/ =
                     DOpEtypes::RefinementType::global)
        {
          //assert(ref_type == DOpEtypes::RefinementType::global);
          RefinementContainer ref_con_dummy;
          RefineTime(ref_con_dummy);
        }

        /******************************************************/
        void
        RefineTime(const RefinementContainer& ref_container)
        {
          DOpEtypes::RefinementType ref_type = ref_container.GetRefType();

          //make sure that we do not use any coarsening
          assert(!ref_container.UsesCoarsening());
          assert(time_triangulation_ != NULL);

          if (DOpEtypes::RefinementType::global == ref_type)
          {
            time_triangulation_->set_all_refine_flags();
          }
          else
          {
            throw DOpEException("Not implemented for name =" + DOpEtypesToString(ref_type),
                "MethodOfLines_SpaceTimeHandler::RefineTime");
          }
          time_triangulation_->prepare_coarsening_and_refinement();

          time_triangulation_->execute_coarsening_and_refinement();
          ReInitTime();
        }
        /******************************************************/

  protected:
    void IncrementStateTicket()
    {
      assert( state_ticket_ < std::numeric_limits<unsigned int>::max());
      state_ticket_++;
    }

    void IncrementControlTicket()
    {
      assert( control_ticket_ < std::numeric_limits<unsigned int>::max());
      control_ticket_++;
    }

  private:
    mutable TimeDoFHandler tdfh_;//FIXME Is it really necessary for tdfh_ and interval_ to be mutable? this is really ugly
    mutable TimeIterator interval_;
    dealii::Triangulation<1>* time_triangulation_;
    unsigned int control_ticket_;
    unsigned int state_ticket_;
    mutable DOpEtypes::ControlType control_type_;
};

}

#endif