spacetimehandler.h

Go to the documentation of this file.

#ifndef SPACE_TIME_HANDLER_H_
#define SPACE_TIME_HANDLER_H_

#include "spacetimehandler_base.h"
#include "active_fe_index_setter_interface.h"
#include "mapping_wrapper.h"
#include "dataout_wrapper.h"

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


// Multi-level routines 
#include <deal.II/multigrid/mg_constrained_dofs.h>
#include <deal.II/multigrid/multigrid.h>
#include <deal.II/multigrid/mg_transfer.h>
#include <deal.II/multigrid/mg_tools.h>
#include <deal.II/multigrid/mg_coarse.h>
#include <deal.II/multigrid/mg_smoother.h>
#include <deal.II/multigrid/mg_matrix.h>



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

namespace DOpE
{
  template<template<int, int> class FE, template<int, int> class DH, typename SPARSITYPATTERN,
      typename VECTOR, int dopedim, int dealdim>
    class SpaceTimeHandler : public SpaceTimeHandlerBase<VECTOR>
    {
      public:
        SpaceTimeHandler(DOpEtypes::ControlType type) :
            SpaceTimeHandlerBase<VECTOR>(type), control_index_(
                dealii::numbers::invalid_unsigned_int), state_index_(
                dealii::numbers::invalid_unsigned_int)
        {
        }
        SpaceTimeHandler(dealii::Triangulation<1> & times,
            DOpEtypes::ControlType type) :
            SpaceTimeHandlerBase<VECTOR>(times, type), control_index_(
                dealii::numbers::invalid_unsigned_int), state_index_(
                dealii::numbers::invalid_unsigned_int)
        {
        }
        SpaceTimeHandler(DOpEtypes::ControlType type,
            const ActiveFEIndexSetterInterface<dopedim, dealdim>& index_setter) :
            SpaceTimeHandlerBase<VECTOR>(type), control_index_(
                dealii::numbers::invalid_unsigned_int), state_index_(
                dealii::numbers::invalid_unsigned_int), fe_index_setter_(
                &index_setter)
        {
        }
        SpaceTimeHandler(dealii::Triangulation<1> & times,
            DOpEtypes::ControlType type,
            const ActiveFEIndexSetterInterface<dopedim, dealdim>& index_setter) :
            SpaceTimeHandlerBase<VECTOR>(times, type), control_index_(
                dealii::numbers::invalid_unsigned_int), state_index_(
                dealii::numbers::invalid_unsigned_int), fe_index_setter_(
                &index_setter)
        {
        }
        virtual
        ~SpaceTimeHandler()
        {

        }

        virtual void
        ReInit(unsigned int control_n_blocks,
            const std::vector<unsigned int>& control_block_component,
            unsigned int state_n_blocks,
            const std::vector<unsigned int>& state_block_component) =0;

        /******************************************************/

        virtual const DOpEWrapper::DoFHandler<dopedim, DH>&
        GetControlDoFHandler() const =0;

        /******************************************************/

        virtual const DOpEWrapper::DoFHandler<dealdim, DH>&
        GetStateDoFHandler() const = 0;

        /******************************************************/

        virtual const DOpEWrapper::Mapping<dealdim, DH>&
        GetMapping() const = 0;

        /******************************************************/

        const std::vector<const DOpEWrapper::DoFHandler<dealdim, DH>*>&
        GetDoFHandler() const
        {
          assert(state_index_ != dealii::numbers::invalid_unsigned_int);
#if dope_dimension > 0
          assert(control_index_ != dealii::numbers::invalid_unsigned_int);
          domain_dofhandler_vector_[control_index_] = &GetControlDoFHandler();
          domain_dofhandler_vector_[state_index_] = &GetStateDoFHandler();
#else
          domain_dofhandler_vector_[state_index_] = &GetStateDoFHandler();
#endif
          return domain_dofhandler_vector_;
        }

        /******************************************************/

        std::vector<
            typename DOpEWrapper::DoFHandler<dealdim, DH>::active_cell_iterator>
        GetDoFHandlerBeginActive() const
        {
          std::vector<
              typename DOpEWrapper::DoFHandler<dealdim, DH>::active_cell_iterator> ret(
              this->GetDoFHandler().size());
          for (unsigned int dh = 0; dh < this->GetDoFHandler().size(); dh++)
          {
            ret[dh] = this->GetDoFHandler()[dh]->begin_active();
          }
          return ret;
        }

        /******************************************************/

        std::vector<
            typename DOpEWrapper::DoFHandler<dealdim, DH>::active_cell_iterator>
        GetDoFHandlerEnd() const
        {
          std::vector<
              typename DOpEWrapper::DoFHandler<dealdim, DH>::active_cell_iterator> ret(
              this->GetDoFHandler().size());
          for (unsigned int dh = 0; dh < this->GetDoFHandler().size(); dh++)
          {
            ret[dh] = this->GetDoFHandler()[dh]->end();
          }
          return ret;
        }


        /******************************************************/

        std::vector<
            typename DOpEWrapper::DoFHandler<dealdim, DH>::cell_iterator>
        GetDoFHandlerBeginActiveAllLevels() const
        {
          std::vector<
              typename DOpEWrapper::DoFHandler<dealdim, DH>::cell_iterator> ret(
              this->GetDoFHandler().size());
          for (unsigned int dh = 0; dh < this->GetDoFHandler().size(); dh++)
          {
            ret[dh] = this->GetDoFHandler()[dh]->begin_active();
          }
          return ret;
        }

        /******************************************************/

        std::vector<
            typename DOpEWrapper::DoFHandler<dealdim, DH>::cell_iterator>
        GetDoFHandlerEndAllLevels() const
        {
          std::vector<
              typename DOpEWrapper::DoFHandler<dealdim, DH>::cell_iterator> ret(
              this->GetDoFHandler().size());
          for (unsigned int dh = 0; dh < this->GetDoFHandler().size(); dh++)
          {
            ret[dh] = this->GetDoFHandler()[dh]->end();
          }
          return ret;
        }



        /******************************************************/

        void
        SetDoFHandlerOrdering(unsigned int control_index,
            unsigned int state_index)
        {
          control_index_ = control_index;
          state_index_ = state_index;
#if dope_dimension > 0
          {
            assert(( control_index_ ==0 && state_index_ ==1 )||( control_index_ ==1 && state_index_ ==0 ));
            domain_dofhandler_vector_.clear();
            if(domain_dofhandler_vector_.size() != 2)
            {
              domain_dofhandler_vector_.resize(2,NULL);
            }
          }
#else
          {
            assert(state_index_ == 0);
            domain_dofhandler_vector_.clear();
            if (domain_dofhandler_vector_.size() != 1)
            {
              domain_dofhandler_vector_.resize(1, NULL);
            }
          }
#endif
        }
        /******************************************************/

        unsigned int
        GetStateIndex()
        {
          return state_index_;
        }

        /******************************************************/
        const ActiveFEIndexSetterInterface<dopedim, dealdim>&
        GetFEIndexSetter() const
        {
          //makes only sense in the hp case.
          return *fe_index_setter_;
        }

        /******************************************************/
        /*
         * This function sets for every element the right fe index for the state variable.
         * This is only useful in the hp case!
         *
         * @param dof_handler   The dof_handler for which the fe indices have to be set.
         */
        void
        SetActiveFEIndicesState(
            DOpEWrapper::DoFHandler<dealdim, DH>& dof_handler)
        {
          if (dof_handler.NeedIndexSetter()) //with this we distinguish between hp and classic
          {
            for (typename DH<dealdim, dealdim>::active_cell_iterator element =
                dof_handler.begin_active(); element != dof_handler.end(); ++element)
            {
              this->GetFEIndexSetter().SetActiveFEIndexState(element);
            }
          }
        }

        /******************************************************/
        /*
         * This function sets for every element the right fe index for the state variable.
         * This is only useful in the hp case!
         *
         * @param dof_handler   The dof_handler for which the fe indices have to be set.
         */
        void
        SetActiveFEIndicesControl(
            DOpEWrapper::DoFHandler<dopedim, DH>& dof_handler)
        {
          if (dof_handler.NeedIndexSetter())
          {
            for (typename DH<dopedim, dopedim>::active_cell_iterator element =
                dof_handler.begin_active(); element != dof_handler.end(); ++element)
            {
              this->GetFEIndexSetter().SetActiveFEIndexState(element);
            }
          }
        }

        /******************************************************/

        virtual unsigned int
        GetControlDoFsPerBlock(unsigned int b, int time_point = -1) const=0;

        /******************************************************/

        virtual unsigned int
        GetStateDoFsPerBlock(unsigned int b, int time_point = -1) const =0;

        /******************************************************/

        virtual unsigned int
        GetConstraintDoFsPerBlock(std::string name, unsigned int b) const=0;

        /******************************************************/

        virtual const std::vector<unsigned int>&
        GetControlDoFsPerBlock(int time_point = -1) const =0;

        /******************************************************/

        virtual const std::vector<unsigned int>&
        GetStateDoFsPerBlock(int time_point = -1) const =0;

        /******************************************************/
        virtual const std::vector<unsigned int>&
        GetConstraintDoFsPerBlock(std::string name) const = 0;

        /******************************************************/

        virtual const dealii::ConstraintMatrix
        &
        GetControlDoFConstraints() const=0;

        /******************************************************/

        virtual const dealii::ConstraintMatrix
        &
        GetStateDoFConstraints() const=0;

        /*******************************************************/

        virtual const std::vector<dealii::Point<dealdim> >
        &
        GetMapDoFToSupportPoints()=0;

        /******************************************************/

        virtual void
        ComputeControlSparsityPattern(SPARSITYPATTERN & sparsity) const=0;

        /******************************************************/

        virtual void
        ComputeStateSparsityPattern(SPARSITYPATTERN & sparsity) const=0;

        /******************************************************/

        virtual void
          ComputeMGStateSparsityPattern(dealii::MGLevelObject<dealii::BlockSparsityPattern> & /*mg_sparsity_patterns*/,
                                        unsigned int /*n_levels*/) const
        {
           throw DOpEException(
                "Not used for normal DofHandler",
                "StateSpaceTimeHandler.h");
        }

        /******************************************************/

        virtual void
          ComputeMGStateSparsityPattern(dealii::MGLevelObject<dealii::SparsityPattern> & /*mg_sparsity_patterns*/,
                                        unsigned int /*n_levels*/) const
        {
           throw DOpEException(
                "Not used for normal DofHandler",
                "StateSpaceTimeHandler.h");
        }

        /******************************************************/
        virtual const FE<dealdim, dealdim>&
        GetFESystem(std::string name) const=0;

        /******************************************************/

        DOpEWrapper::DataOut<dealdim, DH>&
        GetDataOut()
        {
          data_out_.clear();
          return data_out_;
        }

        /******************************************************/

      protected:
        //we need this here, because we know the type of the DoFHandler in use.
        //This saves us a template argument for statpdeproblem etc.
        DOpEWrapper::DataOut<dealdim, DH> data_out_;
        unsigned int control_index_, state_index_;
        const ActiveFEIndexSetterInterface<dopedim, dealdim>* fe_index_setter_;
        mutable std::vector<const DOpEWrapper::DoFHandler<dealdim, DH>*> domain_dofhandler_vector_;
        //TODO What if control and state have different dofhandlertypes??

    };
}

#endif