mol_multimesh_spacetimehandler.h

Go to the documentation of this file.

#ifndef MOL_MULTIMESH_SPACE_TIME_HANDLER_H_
#define MOL_MULTIMESH_SPACE_TIME_HANDLER_H_

#include "mol_spacetimehandler.h"
#include "constraints.h"
#include "sparsitymaker.h"
#include "userdefineddofconstraints.h"
#include "sth_internals.h"
#include "refinementcontainer.h"
#include "solutiontransfer_wrapper.h"

#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/hp/mapping_collection.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/grid/grid_refinement.h>

namespace DOpE
{
  template<template<int, int> class FE, template<int, int> class DH,
      typename SPARSITYPATTERN, typename VECTOR, int dim>
    class MethodOfLines_MultiMesh_SpaceTimeHandler : public SpaceTimeHandler<FE,
        DH, SPARSITYPATTERN, VECTOR, dim, dim>
    {
      public:
        MethodOfLines_MultiMesh_SpaceTimeHandler(
            dealii::Triangulation<dim>& triangulation,
            const FE<dim, dim>& control_fe, const FE<dim, dim>& state_fe,
            DOpEtypes::ControlType type,
            bool flux_pattern = false,
            const ActiveFEIndexSetterInterface<dim, dim>& index_setter =
                ActiveFEIndexSetterInterface<dim, dim>()) :
            SpaceTimeHandler<FE, DH, SPARSITYPATTERN, VECTOR, dim, dim>(type,
                index_setter), state_triangulation_(triangulation), control_dof_handler_(
                control_triangulation_), state_dof_handler_(
                state_triangulation_), control_fe_(&control_fe), state_fe_(
                &state_fe), mapping_(
                DOpEWrapper::StaticMappingQ1<dim, DH>::mapping_q1), constraints_(), control_mesh_transfer_(
                NULL), state_mesh_transfer_(NULL), sparse_mkr_dynamic_(true)
        {
          control_triangulation_.copy_triangulation(state_triangulation_);
          sparsitymaker_ = new SparsityMaker<DH, dim>(flux_pattern);
          user_defined_dof_constr_ = NULL;
        }

        MethodOfLines_MultiMesh_SpaceTimeHandler(
            dealii::Triangulation<dim>& triangulation,
            const FE<dim, dim>& control_fe, const FE<dim, dim>& state_fe,
            dealii::Triangulation<1> & times, DOpEtypes::ControlType type,
            bool flux_pattern = false,
            const ActiveFEIndexSetterInterface<dim, dim>& index_setter =
                ActiveFEIndexSetterInterface<dim, dim>()) :
            SpaceTimeHandler<FE, DH, SPARSITYPATTERN, VECTOR, dim, dim>(times,
                type, index_setter), state_triangulation_(triangulation), control_dof_handler_(
                control_triangulation_), state_dof_handler_(
                state_triangulation_), control_fe_(&control_fe), state_fe_(
                &state_fe), mapping_(
                DOpEWrapper::StaticMappingQ1<dim, DH>::mapping_q1), constraints_(), control_mesh_transfer_(
                NULL), state_mesh_transfer_(NULL), sparse_mkr_dynamic_(true)
        {
          control_triangulation_.copy_triangulation(state_triangulation_);
          sparsitymaker_ = new SparsityMaker<DH, dim>(flux_pattern);
          user_defined_dof_constr_ = NULL;
        }

        MethodOfLines_MultiMesh_SpaceTimeHandler(
            dealii::Triangulation<dim>& triangulation,
            const FE<dim, dim>& control_fe, const FE<dim, dim>& state_fe,
            const Constraints& c, DOpEtypes::ControlType type,
            bool flux_pattern = false,
            const ActiveFEIndexSetterInterface<dim, dim>& index_setter =
                ActiveFEIndexSetterInterface<dim, dim>()) :
            SpaceTimeHandler<FE, DH, SPARSITYPATTERN, VECTOR, dim, dim>(type,
                index_setter), state_triangulation_(triangulation), control_dof_handler_(
                control_triangulation_), state_dof_handler_(
                state_triangulation_), control_fe_(&control_fe), state_fe_(
                &state_fe), mapping_(
                DOpEWrapper::StaticMappingQ1<dim, DH>::mapping_q1), constraints_(
                c), control_mesh_transfer_(NULL), state_mesh_transfer_(NULL), sparse_mkr_dynamic_(
                true)
        {
          control_triangulation_.copy_triangulation(state_triangulation_);
          sparsitymaker_ = new SparsityMaker<DH, dim>(flux_pattern);
          user_defined_dof_constr_ = NULL;
        }

        MethodOfLines_MultiMesh_SpaceTimeHandler(
            dealii::Triangulation<dim>& triangulation,
            const FE<dim, dim>& control_fe, const FE<dim, dim>& state_fe,
            dealii::Triangulation<1> & times, const Constraints& c,
            DOpEtypes::ControlType type,
            bool flux_pattern = false,
            const ActiveFEIndexSetterInterface<dim, dim>& index_setter =
                ActiveFEIndexSetterInterface<dim, dim>()) :
            SpaceTimeHandler<FE, DH, SPARSITYPATTERN, VECTOR, dim, dim>(times,
                type, index_setter), state_triangulation_(triangulation), control_dof_handler_(
                control_triangulation_), state_dof_handler_(
                state_triangulation_), control_fe_(&control_fe), state_fe_(
                &state_fe), mapping_(
                DOpEWrapper::StaticMappingQ1<dim, DH>::mapping_q1), constraints_(
                c), control_mesh_transfer_(NULL), state_mesh_transfer_(NULL), sparse_mkr_dynamic_(
                true)
        {
          control_triangulation_.copy_triangulation(state_triangulation_);
          sparsitymaker_ = new SparsityMaker<DH, dim>(flux_pattern);
          user_defined_dof_constr_ = NULL;
        }

        ~MethodOfLines_MultiMesh_SpaceTimeHandler()
        {
          control_dof_handler_.clear();

          state_dof_handler_.clear();

          if (control_mesh_transfer_ != NULL)
          {
            delete control_mesh_transfer_;
          }
          if (state_mesh_transfer_ != NULL)
          {
            delete state_mesh_transfer_;
          }
          if (sparsitymaker_ != NULL && sparse_mkr_dynamic_ == true)
          {
            delete sparsitymaker_;
          }
        }

        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)
        {
          SpaceTimeHandler<FE, DH, SPARSITYPATTERN, VECTOR, dim, dim>::SetActiveFEIndicesControl(
              control_dof_handler_);
          control_dof_handler_.distribute_dofs(*control_fe_);
          DoFRenumbering::component_wise(static_cast<DH<dim, dim>&>(control_dof_handler_));

          control_dof_constraints_.clear();
          DoFTools::make_hanging_node_constraints(
              static_cast<DH<dim, dim>&>(control_dof_handler_), control_dof_constraints_);
          if (GetUserDefinedDoFConstraints() != NULL)
            GetUserDefinedDoFConstraints()->MakeControlDoFConstraints(
                control_dof_handler_, control_dof_constraints_);
          control_dof_constraints_.close();

          control_dofs_per_block_.resize(control_n_blocks);
          {
            DoFTools::count_dofs_per_block(
                static_cast<DH<dim, dim>&>(control_dof_handler_), control_dofs_per_block_,
                control_block_component);
          }

          SpaceTimeHandler<FE, DH, SPARSITYPATTERN, VECTOR, dim, dim>::SetActiveFEIndicesState(
              state_dof_handler_);
          state_dof_handler_.distribute_dofs(GetFESystem("state"));
          DoFRenumbering::component_wise(static_cast<DH<dim, dim>&>(state_dof_handler_));

          state_dof_constraints_.clear();
          DoFTools::make_hanging_node_constraints(
              static_cast<DH<dim, dim>&>(state_dof_handler_), state_dof_constraints_);
          //TODO Dirichlet Daten hierueber.
          if (GetUserDefinedDoFConstraints() != NULL)
            GetUserDefinedDoFConstraints()->MakeStateDoFConstraints(
                state_dof_handler_, state_dof_constraints_);
          state_dof_constraints_.close();

          state_dofs_per_block_.resize(state_n_blocks);
          DoFTools::count_dofs_per_block(static_cast<DH<dim, dim>&>(state_dof_handler_),
              state_dofs_per_block_, state_block_component);

          support_points_.clear();

          constraints_.ReInit(control_dofs_per_block_);
          //constraints_.ReInit(control_dofs_per_block_, state_dofs_per_block_);

          //Initialize also the timediscretization.
          this->ReInitTime();

          //There where changes invalidate tickets
          this->IncrementControlTicket();
          this->IncrementStateTicket();
        }

        const DOpEWrapper::DoFHandler<dim, DH>&
        GetControlDoFHandler() const
        {
          //There is only one mesh, hence always return this
          return control_dof_handler_;
        }
        const DOpEWrapper::DoFHandler<dim, DH>&
        GetStateDoFHandler() const
        {
          //There is only one mesh, hence always return this
          return state_dof_handler_;
        }

        const DOpEWrapper::Mapping<dim, DH>&
        GetMapping() const
        {
          return mapping_;
        }

        unsigned int
        GetControlDoFsPerBlock(unsigned int b, int /*time_point*/= -1) const
        {
          return control_dofs_per_block_[b];
        }
        unsigned int
        GetStateDoFsPerBlock(unsigned int b, int /*time_point*/= -1) const
        {
          return state_dofs_per_block_[b];
        }
        unsigned int
        GetConstraintDoFsPerBlock(std::string name, unsigned int b) const
        {
          return (constraints_.GetDoFsPerBlock(name))[b];
        }
        const std::vector<unsigned int>&
        GetControlDoFsPerBlock(int /*time_point*/= -1) const
        {
          return control_dofs_per_block_;
        }
        const std::vector<unsigned int>&
        GetStateDoFsPerBlock(int /*time_point*/= -1) const
        {
          return state_dofs_per_block_;
        }
        const std::vector<unsigned int>&
        GetConstraintDoFsPerBlock(std::string name) const
        {
          return constraints_.GetDoFsPerBlock(name);
        }
        const dealii::ConstraintMatrix&
        GetControlDoFConstraints() const
        {
          return control_dof_constraints_;
        }
        const dealii::ConstraintMatrix&
        GetStateDoFConstraints() const
        {
          return state_dof_constraints_;
        }

        virtual void
        InterpolateControl(VECTOR& result,
            const std::vector<VECTOR*> & local_vectors, double t,
            const TimeIterator& it) const
        {
          assert(it.get_left() <= t);
          assert(it.get_right() >= t);
          if (local_vectors.size() != 2)
            throw DOpEException(
                "This function is currently not implemented for anything other than"
                    " linear interpolation of 2 DoFs.",
                "MethodOfLine_SpaceTimeHandler::InterpolateControl");

          double lambda_l = (it.get_right() - t) / it.get_k();
          double lambda_r = (t - it.get_left()) / it.get_k();

          //Here we assume that the numbering of dofs goes from left to right!
          result = *local_vectors[0];

          result.sadd(lambda_l, lambda_r, *local_vectors[1]);
        }
        
        virtual void
        InterpolateState(VECTOR& result,
            const std::vector<VECTOR*> & local_vectors, double t,
            const TimeIterator& it) const
        {
          assert(it.get_left() <= t);
          assert(it.get_right() >= t);
          if (local_vectors.size() != 2)
            throw DOpEException(
                "This function is currently not implemented for anything other than"
                    " linear interpolation of 2 DoFs.",
                "MethodOfLine_SpaceTimeHandler::InterpolateState");

          double lambda_l = (it.get_right() - t) / it.get_k();
          double lambda_r = (t - it.get_left()) / it.get_k();

          //Here we assume that the numbering of dofs goes from left to right!
          result = *local_vectors[0];

          result.sadd(lambda_l, lambda_r, *local_vectors[1]);
        }

        unsigned int
        GetControlNDoFs(int /*time_point*/= -1) const
        {
          return GetControlDoFHandler().n_dofs();
        }
        unsigned int
        GetStateNDoFs(int /*time_point*/= -1) const
        {
          return GetStateDoFHandler().n_dofs();
        }
        unsigned int
        GetConstraintNDoFs(std::string name) const
        {
          return constraints_.n_dofs(name);
        }
        unsigned int
        GetNGlobalConstraints() const
        {
          return constraints_.n_dofs("global");
          //return constraints_.global_n_dofs();
        }
        unsigned int
        GetNLocalConstraints() const
        {
          return constraints_.n_dofs("local");
          //return constraints_.local_n_dofs();
        }

        const std::vector<Point<dim> >&
        GetMapDoFToSupportPoints()
        {
          support_points_.resize(GetStateNDoFs());
          DOpE::STHInternals::MapDoFsToSupportPoints(this->GetMapping(),
              GetStateDoFHandler(), support_points_);
          return support_points_;
        }

        /******************************************************/
        void
        ComputeControlSparsityPattern(SPARSITYPATTERN & sparsity) const;

        /******************************************************/
        void
        ComputeStateSparsityPattern(SPARSITYPATTERN & sparsity) const
        {
          this->GetSparsityMaker()->ComputeSparsityPattern(
              this->GetStateDoFHandler(), sparsity,
              this->GetStateDoFConstraints(), this->GetStateDoFsPerBlock());
        }

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

        const FE<dim, dim>&
        GetFESystem(std::string name) const
        {
          if (name == "state")
          {
            return *state_fe_;
          }
          else if (name == "control")
          {
            return *control_fe_;
          }
          else
          {
            throw DOpEException("Not implemented for name =" + name,
                "MethodOfLines_MultiMesh_SpaceTimeHandler::GetFESystem");
          }

        }

        void
        RefineSpace(DOpEtypes::RefinementType ref_type =
            DOpEtypes::RefinementType::global)
        {
          assert(ref_type == DOpEtypes::RefinementType::global);
          RefinementContainer ref_con_dummy;
          RefineStateSpace(ref_con_dummy);
          RefineControlSpace(ref_con_dummy);
        }

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

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

        template<typename NUMBER>
          void
          RefineSpace(const RefinementContainer&ref_container)
          {
            RefineStateSpace(ref_container);
            RefineControlSpace(ref_container);
          }

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

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

          if (state_mesh_transfer_ != NULL)
          {
            delete state_mesh_transfer_;
            state_mesh_transfer_ = NULL;
          }
          state_mesh_transfer_ = new DOpEWrapper::SolutionTransfer<dim, VECTOR,
            DH>(state_dof_handler_);

          if (DOpEtypes::RefinementType::global == ref_type)
          {
            state_triangulation_.set_all_refine_flags();
          }
          else if (DOpEtypes::RefinementType::fixed_number == ref_type)
          {

            GridRefinement::refine_and_coarsen_fixed_number(
                state_triangulation_, ref_container.GetLocalErrorIndicators(),
                ref_container.GetTopFraction(),
                ref_container.GetBottomFraction());
          }
          else if (DOpEtypes::RefinementType::fixed_fraction == ref_type)
          {

            GridRefinement::refine_and_coarsen_fixed_fraction(
                state_triangulation_, ref_container.GetLocalErrorIndicators(),
                ref_container.GetTopFraction(),
                ref_container.GetBottomFraction());
          }
          else if (DOpEtypes::RefinementType::optimized == ref_type)
          {

            GridRefinement::refine_and_coarsen_optimize(state_triangulation_,
                ref_container.GetLocalErrorIndicators(),
                ref_container.GetConvergenceOrder());
          }
          else if (DOpEtypes::RefinementType::finest_of_both == ref_type)
          {
            this->FlagIfLeftIsNotFinest(state_triangulation_,
                control_triangulation_);
          }
          else
          {
            throw DOpEException("Not implemented for name =" + DOpEtypesToString(ref_type),
                "MethodOfLines_MultiMesh_SpaceTimeHandler::RefineSpace");
          }
          state_triangulation_.prepare_coarsening_and_refinement();
          if (state_mesh_transfer_ != NULL)
            state_mesh_transfer_->prepare_for_pure_refinement();

          state_triangulation_.execute_coarsening_and_refinement();
        }

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

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

          if (control_mesh_transfer_ != NULL)
          {
            delete control_mesh_transfer_;
            control_mesh_transfer_ = NULL;
          }
#if dope_dimension == deal_II_dimension
          control_mesh_transfer_ = new DOpEWrapper::SolutionTransfer<dim, VECTOR,
            DH>(control_dof_handler_);
#endif
          if (DOpEtypes::RefinementType::global == ref_type)
          {
            control_triangulation_.set_all_refine_flags();
          }

          else if (DOpEtypes::RefinementType::fixed_number == ref_type)
          {

            GridRefinement::refine_and_coarsen_fixed_number(
                control_triangulation_, ref_container.GetLocalErrorIndicators(),
                ref_container.GetTopFraction(),
                ref_container.GetBottomFraction());
          }
          else if (DOpEtypes::RefinementType::fixed_fraction == ref_type)
          {

            GridRefinement::refine_and_coarsen_fixed_fraction(
                control_triangulation_, ref_container.GetLocalErrorIndicators(),
                ref_container.GetTopFraction(),
                ref_container.GetBottomFraction());
          }
          else if (DOpEtypes::RefinementType::optimized == ref_type)
          {

            GridRefinement::refine_and_coarsen_optimize(control_triangulation_,
                ref_container.GetLocalErrorIndicators(),
                ref_container.GetConvergenceOrder());
          }
          else if ("finest-of-both")
          {
            this->FlagIfLeftIsNotFinest(control_triangulation_,
                state_triangulation_);
          }
          else
          {
            throw DOpEException("Not implemented for name =" + DOpEtypesToString(ref_type),
                "MethodOfLines_MultiMesh_SpaceTimeHandler::RefineSpace");
          }
          control_triangulation_.prepare_coarsening_and_refinement();

          //FIXME: works only if no coarsening happens, because we do not have the vectors to be interpolated availiable...
          if (control_mesh_transfer_ != NULL)
            control_mesh_transfer_->prepare_for_pure_refinement();

          control_triangulation_.execute_coarsening_and_refinement();
        }
        /******************************************************/

        unsigned int
        NewTimePointToOldTimePoint(unsigned int t) const
        {
          //TODO this has to be implemented when temporal refinement is possible!
          //At present the temporal grid can't be refined
          return t;
        }

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

        void
        SpatialMeshTransferControl(const VECTOR& old_values,
            VECTOR& new_values) const
        {
          if (control_mesh_transfer_ != NULL)
            control_mesh_transfer_->refine_interpolate(old_values, new_values);
        }
        /******************************************************/

        void
        SpatialMeshTransferState(const VECTOR& old_values,
            VECTOR& new_values) const
        {
          if (state_mesh_transfer_ != NULL)
            state_mesh_transfer_->refine_interpolate(old_values, new_values);
        }
        /******************************************************/
        void
        SetUserDefinedDoFConstraints(
            UserDefinedDoFConstraints<DH, dim>& constraints_maker)
        {
          user_defined_dof_constr_ = &constraints_maker;
          user_defined_dof_constr_->RegisterMapping(this->GetMapping());
        }
        /******************************************************/
        void
        SetSparsityMaker(SparsityMaker<DH, dim>& sparsity_maker)
        {
          assert(sparse_mkr_dynamic_==true); //If not true, we already set the sparsity maker
          if (sparsitymaker_ != NULL && sparse_mkr_dynamic_)
            delete sparsitymaker_;
          sparsitymaker_ = &sparsity_maker;
          sparse_mkr_dynamic_ = false;
        }
        /******************************************************/
        void
        ResetStateTriangulation(const dealii::Triangulation<dim>& tria)
        {
          state_dof_handler_.clear();
          state_triangulation_.clear();
          state_triangulation_.copy_triangulation(tria);
          state_dof_handler_.initialize(state_triangulation_, *state_fe_);
          this->IncrementStateTicket();
          if (state_mesh_transfer_ != NULL)
            delete state_mesh_transfer_;
          state_mesh_transfer_ = NULL;
        }

        /******************************************************/
        void
        ResetControlTriangulation(const dealii::Triangulation<dim>& tria)
        {
          control_dof_handler_.clear();
          control_triangulation_.clear();
          control_triangulation_.copy_triangulation(tria);
          control_dof_handler_.initialize(control_triangulation_, *control_fe_);
          this->IncrementControlTicket();
          if (control_mesh_transfer_ != NULL)
            delete control_mesh_transfer_;
          control_mesh_transfer_ = NULL;
        }

      private:
        const SparsityMaker<DH, dim>*
        GetSparsityMaker() const
        {
          return sparsitymaker_;
        }
        const UserDefinedDoFConstraints<DH, dim>*
        GetUserDefinedDoFConstraints() const
        {
          return user_defined_dof_constr_;
        }

        void
        FlagIfLeftIsNotFinest(dealii::Triangulation<dim>& left,
            const dealii::Triangulation<dim>& right);

        SparsityMaker<DH, dim>* sparsitymaker_;
        UserDefinedDoFConstraints<DH, dim>* user_defined_dof_constr_;

        dealii::Triangulation<dim> control_triangulation_;
        dealii::Triangulation<dim>& state_triangulation_;
        DOpEWrapper::DoFHandler<dim, DH> control_dof_handler_;
        DOpEWrapper::DoFHandler<dim, DH> state_dof_handler_;

        std::vector<unsigned int> control_dofs_per_block_;
        std::vector<unsigned int> state_dofs_per_block_;

        dealii::ConstraintMatrix control_dof_constraints_;
        dealii::ConstraintMatrix state_dof_constraints_;

        const dealii::SmartPointer<const FE<dim, dim>> control_fe_;
        const dealii::SmartPointer<const FE<dim, dim>> state_fe_;

        DOpEWrapper::Mapping<dim, DH> mapping_;

        std::vector<Point<dim> > support_points_;

        Constraints constraints_;
        DOpEWrapper::SolutionTransfer<dim, VECTOR,DH>* control_mesh_transfer_;
        DOpEWrapper::SolutionTransfer<dim, VECTOR,DH>* state_mesh_transfer_;
        bool sparse_mkr_dynamic_;
    };

  /**************************explicit instantiation*************/
#if dope_dimension == deal_II_dimension

  template<>
    void
    DOpE::MethodOfLines_MultiMesh_SpaceTimeHandler<dealii::FESystem,
        dealii::DoFHandler, dealii::BlockSparsityPattern,
        dealii::BlockVector<double>, dope_dimension>::ComputeControlSparsityPattern(
        dealii::BlockSparsityPattern & sparsity) const
    {
      const std::vector<unsigned int>& blocks = this->GetControlDoFsPerBlock();
#if DEAL_II_VERSION_GTE(8,3,0)
      dealii::BlockDynamicSparsityPattern csp(blocks.size(),
          blocks.size());
#else
      dealii::BlockCompressedSimpleSparsityPattern csp(blocks.size(),
          blocks.size());
#endif
      for (unsigned int i = 0; i < blocks.size(); i++)
      {
        for (unsigned int j = 0; j < blocks.size(); j++)
        {
          csp.block(i, j).reinit(this->GetControlDoFsPerBlock(i),
              this->GetControlDoFsPerBlock(j));
        }
      }
      csp.collect_sizes();

      dealii::DoFTools::make_sparsity_pattern(
          static_cast<const dealii::DoFHandler<deal_II_dimension>&>(this->GetControlDoFHandler()),
          csp);
      this->GetControlDoFConstraints().condense(csp);
      sparsity.copy_from(csp);
    }

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

  template<>
    void
    MethodOfLines_MultiMesh_SpaceTimeHandler<dealii::FESystem,
        dealii::DoFHandler, dealii::SparsityPattern, dealii::Vector<double>,
        dope_dimension>::ComputeControlSparsityPattern(
        dealii::SparsityPattern & sparsity) const
    {
      const unsigned int total_dofs = this->GetControlNDoFs();
#if DEAL_II_VERSION_GTE(8,3,0)
      dealii::DynamicSparsityPattern csp(total_dofs, total_dofs);
#else
      dealii::CompressedSimpleSparsityPattern csp(total_dofs, total_dofs);
#endif
      dealii::DoFTools::make_sparsity_pattern(
          static_cast<const dealii::DoFHandler<deal_II_dimension>&>(this->GetControlDoFHandler()),
          csp);
      this->GetControlDoFConstraints().condense(csp);
      sparsity.copy_from(csp);
    }

  template<>
    void
    DOpE::MethodOfLines_MultiMesh_SpaceTimeHandler<dealii::hp::FECollection,
        dealii::hp::DoFHandler, dealii::BlockSparsityPattern,
        dealii::BlockVector<double>, dope_dimension>::ComputeControlSparsityPattern(
        dealii::BlockSparsityPattern & sparsity) const
    {
      const std::vector<unsigned int>& blocks = this->GetControlDoFsPerBlock();
#if DEAL_II_VERSION_GTE(8,3,0)
      dealii::BlockDynamicSparsityPattern csp(blocks.size(),
          blocks.size());
#else
      dealii::BlockCompressedSimpleSparsityPattern csp(blocks.size(),
          blocks.size());
#endif      
      for (unsigned int i = 0; i < blocks.size(); i++)
      {
        for (unsigned int j = 0; j < blocks.size(); j++)
        {
          csp.block(i, j).reinit(this->GetControlDoFsPerBlock(i),
              this->GetControlDoFsPerBlock(j));
        }
      }
      csp.collect_sizes();

      dealii::DoFTools::make_sparsity_pattern(
          static_cast<const dealii::hp::DoFHandler<deal_II_dimension>&>(this->GetControlDoFHandler()),
          csp);
      this->GetControlDoFConstraints().condense(csp);
      sparsity.copy_from(csp);
    }

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

  template<>
    void
    MethodOfLines_MultiMesh_SpaceTimeHandler<dealii::hp::FECollection,
        dealii::hp::DoFHandler, dealii::SparsityPattern, dealii::Vector<double>,
        dope_dimension>::ComputeControlSparsityPattern(
        dealii::SparsityPattern & sparsity) const
    {
      const unsigned int total_dofs = this->GetControlNDoFs();
#if DEAL_II_VERSION_GTE(8,3,0)
      dealii::DynamicSparsityPattern csp(total_dofs, total_dofs);
#else
      dealii::CompressedSimpleSparsityPattern csp(total_dofs, total_dofs);
#endif
      dealii::DoFTools::make_sparsity_pattern(
          static_cast<const dealii::hp::DoFHandler<deal_II_dimension>&>(this->GetControlDoFHandler()),
          csp);
      this->GetControlDoFConstraints().condense(csp);
      sparsity.copy_from(csp);
    }
#endif //Endof explicit instanciation
  /*******************************************************/
  template<template<int, int> class FE, template<int, int> class DH,
      typename SPARSITYPATTERN, typename VECTOR, int dim>
    void
    MethodOfLines_MultiMesh_SpaceTimeHandler<FE, DH, SPARSITYPATTERN, VECTOR,
        dim>::FlagIfLeftIsNotFinest(dealii::Triangulation<dim>& left,
        const dealii::Triangulation<dim>& right)
    {
      auto element_list = GridTools::get_finest_common_cells(left, right);
      auto element_iter = element_list.begin();
      for (; element_iter != element_list.end(); element_iter++)
      {
        if (element_iter->second->has_children())
        {
          //left is not finest, so we should flag the left element to be refined
          element_iter->first->set_refine_flag();
        }
      }
    }
}

#endif