integrator.h

Go to the documentation of this file.

#ifndef _Integrator_H_
#define _Integrator_H_

#include <lac/vector.h>
#include <lac/block_sparsity_pattern.h>
#include <lac/block_sparse_matrix.h>

#include <numerics/matrix_tools.h>
#include <numerics/vector_tools.h>

#include <base/function.h>

#include <vector>

#include "elementdatacontainer.h"
#include "facedatacontainer.h"
#include "dwrdatacontainer.h"
#include "residualestimator.h"
#include "dopetypes.h"

namespace DOpE
{
  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    class Integrator
    {
      public:
        Integrator(INTEGRATORDATACONT& idc1);
        Integrator(INTEGRATORDATACONT& idc1, INTEGRATORDATACONT& idc2);

        ~Integrator();

        void
        ReInit();

        template<typename PROBLEM>
          void
          ComputeNonlinearResidual(PROBLEM& pde, VECTOR &residual,
              bool apply_boundary_values = true);
        template<typename PROBLEM>
          void
          ComputeNonlinearLhs(PROBLEM& pde, VECTOR &residual,
              bool apply_boundary_values = true);
        template<typename PROBLEM>
          void
          ComputeNonlinearRhs(PROBLEM& pde, VECTOR &residual,
              bool apply_boundary_values = true);
        template<typename PROBLEM, typename MATRIX>
          void
          ComputeMatrix(PROBLEM& pde, MATRIX &matrix);

         template<typename PROBLEM>
          void
          ComputeNonlinearAlgebraicResidual(PROBLEM& pde, VECTOR &residual);
        template<typename PROBLEM>
          void
          ComputeLocalControlConstraints(PROBLEM& pde, VECTOR &constraints);

        template<typename PROBLEM>
          SCALAR
          ComputeDomainScalar(PROBLEM& pde);
        template<typename PROBLEM>
          SCALAR
          ComputePointScalar(PROBLEM& pde);
        template<typename PROBLEM>
          SCALAR
          ComputeBoundaryScalar(PROBLEM& pde);
        template<typename PROBLEM>
          SCALAR
          ComputeFaceScalar(PROBLEM& pde);
         template<typename PROBLEM>
          SCALAR
          ComputeAlgebraicScalar(PROBLEM& pde);

        template<typename PROBLEM>
          void
          ApplyInitialBoundaryValues(PROBLEM& pde, VECTOR &u);

        template<typename PROBLEM>
          void
          ApplyTransposedInitialBoundaryValues(PROBLEM& pde, VECTOR &u);
        template<typename PROBLEM>
          void
          ApplyNewtonBoundaryValues(PROBLEM& pde, VECTOR &u);
        template<typename PROBLEM, typename MATRIX>
          void
          ApplyNewtonBoundaryValues(PROBLEM& pde, MATRIX &matrix, VECTOR &rhs,
              VECTOR &sol);

        inline void
        AddDomainData(std::string name, const VECTOR* new_data);
        inline void
        DeleteDomainData(std::string name);
        inline void
        AddParamData(std::string name, const dealii::Vector<SCALAR>* new_data);inline void
        DeleteParamData(std::string name);
 
        template<typename PROBLEM, class STH, class EDC, class FDC>
          void
          ComputeRefinementIndicators(PROBLEM& pde,
              DWRDataContainer<STH, INTEGRATORDATACONT, EDC, FDC, VECTOR>& dwrc);
        template<typename PROBLEM>
          void
          ComputeRefinementIndicators(PROBLEM& pde,
              ResidualErrorContainer<VECTOR>& dwrc);

        inline INTEGRATORDATACONT&
        GetIntegratorDataContainer() const;

        inline INTEGRATORDATACONT&
        GetIntegratorDataContainerFunc() const;

      protected: 
        inline const std::map<std::string, const VECTOR*>&
          GetDomainData() const;
        inline const std::map<std::string, const dealii::Vector<SCALAR>*>&
          GetParamData() const;

      private:
        template<template<int, int> class DH>
          void
          InterpolateBoundaryValues(
              const DOpEWrapper::Mapping<dim, DH>& mapping,
              const DOpEWrapper::DoFHandler<dim, DH>* dof_handler,
              const unsigned int color, const dealii::Function<dim>& function,
              std::map<unsigned int, SCALAR>& boundary_values,
              const std::vector<bool>& comp_mask) const;

        template<typename ELEMENTITERATOR>
          bool
          AtInterface(ELEMENTITERATOR& element, unsigned int face)
          {
            if (element[0]->neighbor_index(face) != -1)
              if (element[0]->material_id()
                  != element[0]->neighbor(face)->material_id())
                return true;
            return false;
          }

        INTEGRATORDATACONT & _idc1;
        INTEGRATORDATACONT & _idc2;

        std::map<std::string, const VECTOR*> _domain_data;
        std::map<std::string, const dealii::Vector<SCALAR>*> _param_data;
    };

  /**********************************Implementation*******************************************/

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::Integrator(
        INTEGRATORDATACONT& idc)
        : _idc1(idc), _idc2(idc)
    {
    }

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::Integrator(
        INTEGRATORDATACONT& idc1, INTEGRATORDATACONT& idc2)
        : _idc1(idc1), _idc2(idc2)
    {
    }

  /**********************************Implementation*******************************************/

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::~Integrator()
    {

    }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    void
    Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ReInit()
    {

    }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<typename PROBLEM>
      void
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeNonlinearResidual(
          PROBLEM& pde, VECTOR &residual, bool apply_boundary_values)
      {
        residual = 0.;
        // Begin integration
        unsigned int dofs_per_element;
        dealii::Vector<SCALAR> local_vector;
        std::vector<unsigned int> local_dof_indices;

        const bool need_point_rhs = pde.HasPoints();

        const auto& dof_handler =
            pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
        auto element =
            pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
        auto endc =
            pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();

        // Generate the data containers.
        GetIntegratorDataContainer().InitializeEDC(pde.GetUpdateFlags(),
            *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
            this->GetParamData(), this->GetDomainData());
        auto& edc = GetIntegratorDataContainer().GetElementDataContainer();

        bool need_faces = pde.HasFaces();
        bool need_interfaces = pde.HasInterfaces();
        std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
        bool need_boundary_integrals = (boundary_equation_colors.size() > 0);
        GetIntegratorDataContainer().InitializeFDC(pde.GetFaceUpdateFlags(),
            *(pde.GetBaseProblem().GetSpaceTimeHandler()),
            element,
            this->GetParamData(),
            this->GetDomainData(),
            need_interfaces);
        auto & fdc = GetIntegratorDataContainer().GetFaceDataContainer();

        for (; element[0] != endc[0]; element[0]++)
        {
          for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
          {
            if (element[dh] == endc[dh])
            {
              throw DOpEException(
                  "Elementnumbers in DoFHandlers are not matching!",
                  "Integrator::ComputeNonlinearResidual");
            }
          }

          edc.ReInit();

          dofs_per_element = element[0]->get_fe().dofs_per_cell;

          local_vector.reinit(dofs_per_element);
          local_vector = 0;

          local_dof_indices.resize(0);
          local_dof_indices.resize(dofs_per_element, 0);

          //the second '1' plays only a role in the stationary case. In the non-stationary
          //case, scale_ico is set by the time-stepping-scheme
          pde.ElementEquation(edc, local_vector, 1., 1.);
          pde.ElementRhs(edc, local_vector, -1.);

          if(need_boundary_integrals && element[0]->at_boundary())
          {
            for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
            {
              if (element[0]->face(face)->at_boundary()
                  &&
                  (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
                          element[0]->face(face)->boundary_indicator()) != boundary_equation_colors.end()))
              {
                fdc.ReInit(face);
                pde.BoundaryEquation(fdc,local_vector, 1., 1.);
                pde.BoundaryRhs(fdc,local_vector,-1.);
              }
            }
          }
          if(need_faces)
          {
            for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
            {
              if (element[0]->neighbor_index(face) != -1)
              {
                fdc.ReInit(face);
                pde.FaceEquation(fdc, local_vector, 1., 1.);
                pde.FaceRhs(fdc, local_vector,-1.);
              }
            }
          }
          if( need_interfaces)
          {

            for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
            {
              // auto face_it = element[0]->face(face);
              // first, check if we are at an interface, i.e. not the neighbour exists and
              // it has a different material_id than the actual element
              if(AtInterface(element, face))
              {
                //There exist now 3 different scenarios, given the actual element and face:
                // The neighbour behind this face is [ more | as much | less] refined
                // than/as the actual element. We have to distinguish here only between the case 1
                // and the other two, because these will be distinguished in in the FaceDataContainer.

                if (element[0]->neighbor(face)->has_children())
                {
                  //first: neighbour is finer

                  for (unsigned int subface_no=0;
                      subface_no < element[0]->face(face)->n_children();
                      ++subface_no)
                  {
                    //TODO Now here we have to initialise the subface_values on the
                    // actual element and then the facevalues of the neighbours
                    fdc.ReInit(face, subface_no);
                    fdc.ReInitNbr();

                    pde.InterfaceEquation(fdc, local_vector, 1., 1.);

                  }
                }
                else
                {
                  // either neighbor is as fine as this element or
                  // it is coarser

                  fdc.ReInit(face);
                  fdc.ReInitNbr();
                  pde.InterfaceEquation(fdc, local_vector, 1., 1.);
                }

              }                  //endif atinterface
            }                  //endfor faces
          }                  //endif need_interfaces
          //LocalToGlobal
          element[0]->get_dof_indices(local_dof_indices);
          for (unsigned int i = 0; i < dofs_per_element; ++i)
          {
            residual(local_dof_indices[i]) += local_vector(i);
          }

          for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
          {
            element[dh]++;
          }
        }

        //check if we need the evaluation of PointRhs
        if (need_point_rhs)
        {
          VECTOR point_rhs;
          pde.PointRhs(this->GetParamData(), this->GetDomainData(), point_rhs,
              -1.);
          residual += point_rhs;
        }

        if (apply_boundary_values)
        {
          ApplyNewtonBoundaryValues(pde, residual);
        }
      }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<typename PROBLEM>
      void
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeNonlinearLhs(
          PROBLEM& pde, VECTOR &residual, bool apply_boundary_values)
      {
        {

          residual = 0.;
          // Begin integration
          unsigned int dofs_per_element;

          dealii::Vector<SCALAR> local_vector;

          std::vector<unsigned int> local_dof_indices;

          const auto& dof_handler =
              pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
          auto element =
              pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
          auto endc =
              pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();

          // Generate the data containers.
          GetIntegratorDataContainer().InitializeEDC(pde.GetUpdateFlags(),
              *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
              this->GetParamData(), this->GetDomainData());
          auto& edc = GetIntegratorDataContainer().GetElementDataContainer();

          bool need_faces = pde.HasFaces();
          bool need_interfaces = pde.HasInterfaces();
          std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
          bool need_boundary_integrals = (boundary_equation_colors.size() > 0);

          GetIntegratorDataContainer().InitializeFDC(pde.GetFaceUpdateFlags(),
              *(pde.GetBaseProblem().GetSpaceTimeHandler()),
              element,
              this->GetParamData(),
              this->GetDomainData(),
              need_interfaces);
          auto & fdc = GetIntegratorDataContainer().GetFaceDataContainer();

          for (; element[0] != endc[0]; element[0]++)
          {
            for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
            {
              if (element[dh] == endc[dh])
              {
                throw DOpEException(
                    "Elementnumbers in DoFHandlers are not matching!",
                    "mIntegrator::ComputeNonlinearLhs");
              }
            }

            edc.ReInit();
            dofs_per_element = element[0]->get_fe().dofs_per_cell;

            local_vector.reinit(dofs_per_element);
            local_vector = 0;

            local_dof_indices.resize(0);
            local_dof_indices.resize(dofs_per_element, 0);

            //the second '1' plays only a role in the stationary case. In the non-stationary
            //case, scale_ico is set by the time-stepping-scheme
            pde.ElementEquation(edc, local_vector, 1., 1.);

            if(need_boundary_integrals)
            {
              for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
              {
                if (element[0]->face(face)->at_boundary()
                    &&
                    (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
                            element[0]->face(face)->boundary_indicator()) != boundary_equation_colors.end()))
                {
                  fdc.ReInit(face);
                  pde.BoundaryEquation(fdc,local_vector, 1., 1.);
                }
              }
            }
            if(need_faces)
            {
              for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
              {
                if (element[0]->neighbor_index(face) != -1)
                {
                  fdc.ReInit(face);
                  pde.FaceEquation(fdc, local_vector, 1., 1.);
                }
              }
            }

            if( need_interfaces)
            {

              for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
              {
                //auto face_it = element[0]->face(face);
                // first, check if we are at an interface, i.e. not the neighbour exists and
                // it has a different material_id than the actual element
                if(AtInterface(element, face))
                {
                  //There exist now 3 different scenarios, given the actual element and face:
                  // The neighbour behind this face is [ more | as much | less] refined
                  // than/as the actual element. We have to distinguish here only between the case 1
                  // and the other two, because these will be distinguished in in the FaceDataContainer.

                  if (element[0]->neighbor(face)->has_children())
                  {
                    //first: neighbour is finer

                    for (unsigned int subface_no=0;
                        subface_no < element[0]->face(face)->n_children();
                        ++subface_no)
                    {
                      //TODO Now here we have to initialise the subface_values on the
                      // actual element and then the facevalues of the neighbours
                      fdc.ReInit(face, subface_no);
                      fdc.ReInitNbr();

                      pde.InterfaceEquation(fdc, local_vector, 1., 1.);

                    }
                  }
                  else
                  {
                    // either neighbor is as fine as this element or
                    // it is coarser

                    fdc.ReInit(face);
                    fdc.ReInitNbr();
                    pde.InterfaceEquation(fdc, local_vector, 1., 1.);
                  }
                }                    //endif atinterface
              }                    //endfor face
            }                    //endif need_interfaces
            //LocalToGlobal
            element[0]->get_dof_indices(local_dof_indices);
            for (unsigned int i = 0; i < dofs_per_element; ++i)
            {
              residual(local_dof_indices[i]) += local_vector(i);
            }

            for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
            {
              element[dh]++;
            }
          }

          if (apply_boundary_values)
          {
            ApplyNewtonBoundaryValues(pde, residual);
          }
        }
      }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<typename PROBLEM>
      void
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeNonlinearRhs(
          PROBLEM& pde, VECTOR &residual, bool apply_boundary_values)
      {
        residual = 0.;
        // Begin integration
        unsigned int dofs_per_element;
        dealii::Vector<SCALAR> local_vector;
        std::vector<unsigned int> local_dof_indices;

        const bool need_point_rhs = pde.HasPoints();

        const auto& dof_handler =
            pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
        auto element =
            pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
        auto endc =
            pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();

        // Initialize the data containers.
        GetIntegratorDataContainer().InitializeEDC(pde.GetUpdateFlags(),
            *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
            this->GetParamData(), this->GetDomainData());
        auto& edc = GetIntegratorDataContainer().GetElementDataContainer();

        bool need_interfaces = pde.HasInterfaces();
        if(need_interfaces )
        {
          throw DOpEException(" Faces on multiple meshes not implemented yet!",
              "IntegratorMultiMesh::ComputeNonlinearRhs_Recursive");
        }
        bool need_faces = pde.HasFaces();
        std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
        bool need_boundary_integrals = (boundary_equation_colors.size() > 0);

        GetIntegratorDataContainer().InitializeFDC(pde.GetFaceUpdateFlags(),
            *(pde.GetBaseProblem().GetSpaceTimeHandler()),
            element,
            this->GetParamData(),
            this->GetDomainData());
        auto & fdc = GetIntegratorDataContainer().GetFaceDataContainer();

        for (; element[0] != endc[0]; element[0]++)
        {
          for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
          {
            if (element[dh] == endc[dh])
            {
              throw DOpEException(
                  "Elementnumbers in DoFHandlers are not matching!",
                  "Integrator::ComputeNonlinearRhs");
            }
          }

          edc.ReInit();
          dofs_per_element = element[0]->get_fe().dofs_per_cell;

          local_vector.reinit(dofs_per_element);
          local_vector = 0;

          local_dof_indices.resize(0);
          local_dof_indices.resize(dofs_per_element, 0);
          pde.ElementRhs(edc, local_vector, 1.);

          if(need_boundary_integrals)
          {
            for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
            {
              if (element[0]->face(face)->at_boundary()
                  &&
                  (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
                          element[0]->face(face)->boundary_indicator()) != boundary_equation_colors.end()))
              {
                fdc.ReInit(face);
                pde.BoundaryRhs(fdc,local_vector,1.);
              }
            }
          }
          if(need_faces)
          {
            for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
            {
              if (element[0]->neighbor_index(face) != -1)
              {
                fdc.ReInit(face);
                pde.FaceRhs(fdc, local_vector);
              }
            }
          }
          //LocalToGlobal
          element[0]->get_dof_indices(local_dof_indices);
          for (unsigned int i = 0; i < dofs_per_element; ++i)
          {
            residual(local_dof_indices[i]) += local_vector(i);
          }

          for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
          {
            element[dh]++;
          }
        }

        //check if we need the evaluation of PointRhs
        if (need_point_rhs)
        {
          VECTOR point_rhs;
          pde.PointRhs(this->GetParamData(), this->GetDomainData(), point_rhs,
              1.);
          residual += point_rhs;
        }

        if (apply_boundary_values)
        {
          ApplyNewtonBoundaryValues(pde, residual);
        }
      }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<typename PROBLEM, typename MATRIX>
      void
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeMatrix(
          PROBLEM& pde, MATRIX &matrix)
      {
        matrix = 0.;
        // Begin integration
        unsigned int dofs_per_element;
        std::vector<unsigned int> local_dof_indices;

        const auto& dof_handler =
            pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
        auto element =
            pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
        auto endc =
            pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();

        GetIntegratorDataContainer().InitializeEDC(pde.GetUpdateFlags(),
            *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
            this->GetParamData(), this->GetDomainData());
        auto& edc = GetIntegratorDataContainer().GetElementDataContainer();

        //for the interface-case
        unsigned int nbr_dofs_per_element;
        std::vector<unsigned int> nbr_local_dof_indices;

        bool need_faces = pde.HasFaces();
        bool need_interfaces = pde.HasInterfaces();
        std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
        bool need_boundary_integrals = (boundary_equation_colors.size() > 0);

        GetIntegratorDataContainer().InitializeFDC(pde.GetFaceUpdateFlags(),
            *(pde.GetBaseProblem().GetSpaceTimeHandler()),
            element,
            this->GetParamData(),
            this->GetDomainData(),
            need_interfaces);
        auto & fdc = GetIntegratorDataContainer().GetFaceDataContainer();

        for (; element[0] != endc[0]; element[0]++)
        {
          for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
          {
            if (element[dh] == endc[dh])
            {
              throw DOpEException(
                  "Elementnumbers in DoFHandlers are not matching!",
                  "Integrator::ComputeMatrix");
            }
          }
          edc.ReInit();
          dofs_per_element = element[0]->get_fe().dofs_per_cell;

          dealii::FullMatrix<SCALAR> local_matrix(dofs_per_element,
              dofs_per_element);
          local_matrix = 0;

          local_dof_indices.resize(0);
          local_dof_indices.resize(dofs_per_element, 0);
          pde.ElementMatrix(edc, local_matrix);

          if(need_boundary_integrals)
          {
            for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
            {
              if (element[0]->face(face)->at_boundary()
                  &&
                  (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
                          element[0]->face(face)->boundary_indicator()) != boundary_equation_colors.end()))
              {
                fdc.ReInit(face);
                pde.BoundaryMatrix(fdc, local_matrix);
              }
            }
          }
          if(need_faces)
          {
            for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
            {
              if (element[0]->neighbor_index(face) != -1)
              {
                fdc.ReInit(face);
                pde.FaceMatrix(fdc, local_matrix);
              }
            }
          }

          if( need_interfaces)
          {
            for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
            {
              //auto face_it = element[0]->face(face);
              // first, check if we are at an interface, i.e. not the neighbour exists and
              // it has a different material_id than the actual element
              if(AtInterface(element, face))
              {
                //There exist now 3 different scenarios, given the actual element and face:
                // The neighbour behind this face is [ more | as much | less] refined
                // than/as the actual element. We have to distinguish here only between the case 1
                // and the other two, because these will be distinguished in in the FaceDataContainer.

                if (element[0]->neighbor(face)->has_children())
                {
                  //first: neighbour is finer

                  for (unsigned int subface_no=0;
                      subface_no < element[0]->face(face)->n_children();
                      ++subface_no)
                  {
                    //TODO Now here we have to initialise the subface_values on the
                    // actual element and then the facevalues of the neighbours
                    fdc.ReInit(face, subface_no);
                    fdc.ReInitNbr();

                    //TODO auslagern?
                    nbr_dofs_per_element = fdc.GetNbrNDoFsPerElement();
                    nbr_local_dof_indices.resize(0);
                    nbr_local_dof_indices.resize(nbr_dofs_per_element, 0);
                    dealii::FullMatrix<SCALAR> local_interface_matrix(dofs_per_element,nbr_dofs_per_element );
                    local_interface_matrix = 0;

                    pde.InterfaceMatrix(fdc, local_interface_matrix);

                    element[0]->get_dof_indices(local_dof_indices);
                    element[0]->neighbor(face)->get_dof_indices(nbr_local_dof_indices);

                    for (unsigned int i = 0; i < dofs_per_element; ++i)
                    {
                      for (unsigned int j = 0; j < nbr_dofs_per_element; ++j)
                      {
                        matrix.add(local_dof_indices[i], nbr_local_dof_indices[j],
                            local_interface_matrix(i, j));
                      } //endfor j
                    } //endfor i

                  }
                }
                else
                {
                  // either neighbor is as fine as this element or it is coarser
                  fdc.ReInit(face);
                  fdc.ReInitNbr();

                  //TODO auslagern?
                  nbr_dofs_per_element = fdc.GetNbrNDoFsPerElement();
                  nbr_local_dof_indices.resize(0);
                  nbr_local_dof_indices.resize(nbr_dofs_per_element, 0);
                  dealii::FullMatrix<SCALAR> local_interface_matrix(dofs_per_element,nbr_dofs_per_element );
                  local_interface_matrix = 0;

                  pde.InterfaceMatrix(fdc, local_interface_matrix);

                  element[0]->get_dof_indices(local_dof_indices);
                  element[0]->neighbor(face)->get_dof_indices(nbr_local_dof_indices);

                  for (unsigned int i = 0; i < dofs_per_element; ++i)
                  {
                    for (unsigned int j = 0; j < nbr_dofs_per_element; ++j)
                    {
                      matrix.add(local_dof_indices[i], nbr_local_dof_indices[j],
                          local_interface_matrix(i, j));
                    } //endfor j
                  } //endfor i
                }
              } //endif atinterface

            }
          } //endif need_interfaces

          //LocalToGlobal
          element[0]->get_dof_indices(local_dof_indices);
          for (unsigned int i = 0; i < dofs_per_element; ++i)
          {
            for (unsigned int j = 0; j < dofs_per_element; ++j)
            {
              matrix.add(local_dof_indices[i], local_dof_indices[j],
                  local_matrix(i, j));
            }
          }

          for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
          {
            element[dh]++;
          }
        }
      }

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<typename PROBLEM>
      SCALAR
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeDomainScalar(
          PROBLEM& pde)
      {
        {
          SCALAR ret = 0.;

          const auto& dof_handler =
              pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
          auto element =
              pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
          auto endc =
              pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();
          GetIntegratorDataContainerFunc().InitializeEDC(pde.GetUpdateFlags(),
              *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
              this->GetParamData(), this->GetDomainData());
          auto& edc = GetIntegratorDataContainerFunc().GetElementDataContainer();

          bool need_faces = pde.HasFaces();

          for (; element[0] != endc[0]; element[0]++)
          {
            for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
            {
              if (element[dh] == endc[dh])
              {
                throw DOpEException(
                    "Elementnumbers in DoFHandlers are not matching!",
                    "Integrator::ComputeDomainScalar");
              }
            }

            edc.ReInit();
            ret += pde.ElementFunctional(edc);

            if (need_faces)
            {
              throw DOpEException("Face Integrals not Implemented!",
                  "Integrator::ComputeDomainScalar");
            }

            for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
            {
              element[dh]++;
            }
          }
          return ret;
        }
      }
  /*******************************************************************************************/

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<typename PROBLEM>
      SCALAR
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputePointScalar(
          PROBLEM& pde)
      {

        {
          SCALAR ret = 0.;
          ret += pde.PointFunctional(this->GetParamData(),
              this->GetDomainData());

          return ret;
        }
      }
  /*******************************************************************************************/

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<typename PROBLEM>
      SCALAR
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeBoundaryScalar(
          PROBLEM& pde
          )
      {
        {
          SCALAR ret = 0.;
          // Begin integration
          const auto& dof_handler =
          pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
          auto element = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
          auto endc = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();

          GetIntegratorDataContainerFunc().InitializeFDC(pde.GetFaceUpdateFlags(),
              *(pde.GetBaseProblem().GetSpaceTimeHandler()),
              element,
              this->GetParamData(),
              this->GetDomainData());
          auto & fdc = GetIntegratorDataContainerFunc().GetFaceDataContainer();

          std::vector<unsigned int> boundary_functional_colors = pde.GetBoundaryFunctionalColors();
          bool need_boundary_integrals = (boundary_functional_colors.size() > 0);
          if(!need_boundary_integrals)
          {
            throw DOpEException("No boundary colors given!","Integrator::ComputeBoundaryScalar");
          }

          for (;element[0]!=endc[0]; element[0]++)
          {
            for(unsigned int dh=1; dh<dof_handler.size(); dh++)
            {
              if( element[dh] == endc[dh])
              {
                throw DOpEException("Elementnumbers in DoFHandlers are not matching!","Integrator::ComputeBoundaryScalar");
              }
            }

            if(need_boundary_integrals)
            {
              for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
              {
                if (element[0]->face(face)->at_boundary()
                    &&
                    (find(boundary_functional_colors.begin(),boundary_functional_colors.end(),
                            element[0]->face(face)->boundary_indicator()) != boundary_functional_colors.end()))
                {
                  fdc.ReInit(face);
                  ret += pde.BoundaryFunctional(fdc);
                }
              }
            }
            for(unsigned int dh=1; dh<dof_handler.size(); dh++)
            {
              element[dh]++;
            }
          }

          return ret;

        }
      }
  /*******************************************************************************************/

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<typename PROBLEM>
      SCALAR
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeFaceScalar(
          PROBLEM& pde
          )
      {

        {
          SCALAR ret = 0.;
          // Begin integration
          const auto& dof_handler =
          pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
          auto element = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
          auto endc = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();

          GetIntegratorDataContainerFunc().InitializeFDC(pde.GetFaceUpdateFlags(),
              *(pde.GetBaseProblem().GetSpaceTimeHandler()),
              element,
              this->GetParamData(),
              this->GetDomainData());
          auto & fdc = GetIntegratorDataContainerFunc().GetFaceDataContainer();

          bool need_faces = pde.HasFaces();
          if(!need_faces)
          {
            throw DOpEException("No faces required!","Integrator::ComputeFaceScalar");
          }

          for (;element[0]!=endc[0]; element[0]++)
          {
            for(unsigned int dh=1; dh<dof_handler.size(); dh++)
            {
              if( element[dh] == endc[dh])
              {
                throw DOpEException("Elementnumbers in DoFHandlers are not matching!","Integrator::ComputeFaceScalar");
              }
            }

            if(need_faces)
            {
              for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
              {
                if (element[0]->neighbor_index(face) != -1)
                {
                  fdc.ReInit(face);
                  ret +=pde.FaceFunctional(fdc);
                }
              }
            }
            for(unsigned int dh=1; dh<dof_handler.size(); dh++)
            {
              element[dh]++;
            }
          }
          return ret;
        }
      }
  /*******************************************************************************************/

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<typename PROBLEM>
      SCALAR
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeAlgebraicScalar(
          PROBLEM& pde)
      {

        {
          SCALAR ret = 0.;
          ret = pde.AlgebraicFunctional(this->GetParamData(),
              this->GetDomainData());
          return ret;
        }
      }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<typename PROBLEM>
      void
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeNonlinearAlgebraicResidual(
          PROBLEM& pde, VECTOR &residual)
      {
        residual = 0.;
        pde.AlgebraicResidual(residual, this->GetParamData(),
            this->GetDomainData());
      }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<typename PROBLEM>
      void
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeLocalControlConstraints(
          PROBLEM& pde, VECTOR &constraints)
      {
        constraints = 0.;
        pde.ComputeLocalControlConstraints(constraints, this->GetParamData(),
            this->GetDomainData());
      }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<typename PROBLEM>
      void
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ApplyTransposedInitialBoundaryValues(
        PROBLEM& /*pde*/, VECTOR &/*u*/)
      {
        //Wird  hier nicht gebraucht...
      }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<typename PROBLEM>
      void
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ApplyInitialBoundaryValues(
          PROBLEM& pde, VECTOR &u)
      {
        //TODO Apply constraints locally, see, e.g., dealii step-27 ? But howto do this in the newton iter
        // e.g. sometimes we need zero sometimes we need other values.

        //Never Condense Nodes Here ! Or All will fail if the state is not initialized with zero!
        //pde.GetDoFConstraints().condense(u);
        std::vector<unsigned int> dirichlet_colors = pde.GetDirichletColors();
        for (unsigned int i = 0; i < dirichlet_colors.size(); i++)
        {
          unsigned int color = dirichlet_colors[i];
          std::vector<bool> comp_mask = pde.GetDirichletCompMask(color);
          std::map<unsigned int, SCALAR> boundary_values;

          InterpolateBoundaryValues( pde.GetBaseProblem().GetSpaceTimeHandler()->GetMapping(),
              pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler()[0],
              color,
              pde.GetDirichletValues(color, this->GetParamData(),
                  this->GetDomainData()), boundary_values, comp_mask);

          for (typename std::map<unsigned int, SCALAR>::const_iterator p =
              boundary_values.begin(); p != boundary_values.end(); p++)
          {
            u(p->first) = p->second;
          }
        }
      }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<typename PROBLEM>
      void
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ApplyNewtonBoundaryValues(
          PROBLEM& pde, VECTOR &u)
      {
        //TODO Apply constraints locally, see, e.g., dealii step-27 ? But howto do this in the newton iter
        // e.g. sometimes we need zero sometimes we need other values.

        pde.GetDoFConstraints().condense(u);
        std::vector<unsigned int> dirichlet_colors = pde.GetDirichletColors();
        for (unsigned int i = 0; i < dirichlet_colors.size(); i++)
        {
          unsigned int color = dirichlet_colors[i];
          std::vector<bool> comp_mask = pde.GetDirichletCompMask(color);
          std::map<unsigned int, SCALAR> boundary_values;

          InterpolateBoundaryValues( pde.GetBaseProblem().GetSpaceTimeHandler()->GetMapping(),
              pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler()[0],
              color, dealii::ZeroFunction<dim>(comp_mask.size()),
              boundary_values, comp_mask);

          for (typename std::map<unsigned int, SCALAR>::const_iterator p =
              boundary_values.begin(); p != boundary_values.end(); p++)
          {
            u(p->first) = p->second;
          }
        }
      }
  /*******************************************************************************************/

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<typename PROBLEM, typename MATRIX>
      void
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ApplyNewtonBoundaryValues(
          PROBLEM& pde, MATRIX& matrix, VECTOR &rhs, VECTOR &sol)
      {
        //TODO Apply constraints locally, see, e.g., dealii step-27 ? But howto do this in the newton iter
        // e.g. sometimes we need zero sometimes we need other values.
        pde.GetDoFConstraints().condense(rhs);
        pde.GetDoFConstraints().condense(matrix);
        std::vector<unsigned int> dirichlet_colors = pde.GetDirichletColors();
        for (unsigned int i = 0; i < dirichlet_colors.size(); i++)
        {
          unsigned int color = dirichlet_colors[i];
          std::vector<bool> comp_mask = pde.GetDirichletCompMask(color);
          std::map<unsigned int, SCALAR> boundary_values;

          InterpolateBoundaryValues( pde.GetBaseProblem().GetSpaceTimeHandler()->GetMapping(),
              pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler()[0],
              color, dealii::ZeroFunction<dim>(comp_mask.size()),
              boundary_values, comp_mask);

          dealii::MatrixTools::apply_boundary_values(boundary_values, matrix,
              sol, rhs);
        }
      }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    void
    Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::AddDomainData(
        std::string name, const VECTOR* new_data)
    {
      if (_domain_data.find(name) != _domain_data.end())
      {
        throw DOpEException(
            "Adding multiple Data with name " + name + " is prohibited!",
            "Integrator::AddDomainData");
      }
      _domain_data.insert(
          std::pair<std::string, const VECTOR*>(name, new_data));
    }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    void
    Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::DeleteDomainData(
        std::string name)
    {
      typename std::map<std::string, const VECTOR *>::iterator it =
          _domain_data.find(name);
      if (it == _domain_data.end())
      {
        throw DOpEException(
            "Deleting Data " + name + " is impossible! Data not found",
            "Integrator::DeleteDomainData");
      }
      _domain_data.erase(it);
    }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    const std::map<std::string, const VECTOR*>&
    Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::GetDomainData() const
    {
      return _domain_data;
    }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    void
    Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::AddParamData(
        std::string name, const dealii::Vector<SCALAR>* new_data)
    {
      if (_param_data.find(name) != _param_data.end())
      {
        throw DOpEException(
            "Adding multiple Data with name " + name + " is prohibited!",
            "Integrator::AddParamData");
      }
      _param_data.insert(
          std::pair<std::string, const dealii::Vector<SCALAR>*>(name,
              new_data));
    }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    void
    Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::DeleteParamData(
        std::string name)
    {
      typename std::map<std::string, const dealii::Vector<SCALAR>*>::iterator it =
          _param_data.find(name);
      if (it == _param_data.end())
      {
        throw DOpEException(
            "Deleting Data " + name + " is impossible! Data not found",
            "Integrator::DeleteParamData");
      }
      _param_data.erase(it);
    }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    const std::map<std::string, const dealii::Vector<SCALAR>*>&
    Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::GetParamData() const
    {
      return _param_data;
    }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<typename PROBLEM, class STH, class EDC, class FDC>
      void
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeRefinementIndicators(
          PROBLEM& pde,
          DWRDataContainer<STH, INTEGRATORDATACONT, EDC, FDC, VECTOR>& dwrc)
      {
        unsigned int n_error_comps = dwrc.GetNErrorComps();
        //for primal and dual part of the error
        std::vector<double> element_sum(n_error_comps, 0);
        element_sum.resize(n_error_comps, 0);

        // Begin integration
        const auto& dof_handler =
            pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
        auto element =
            pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
        auto endc =
            pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();

        const auto& dof_handler_weight = dwrc.GetWeightSTH().GetDoFHandler();
        auto element_weight = dwrc.GetWeightSTH().GetDoFHandlerBeginActive();
        auto endc_high = dwrc.GetWeightSTH().GetDoFHandlerEnd();

        // Generate the data containers. Notice that we use the quadrature
        //formula from the higher order idc!.
        GetIntegratorDataContainer().InitializeEDC(
            dwrc.GetWeightIDC().GetQuad(), pde.GetUpdateFlags(),
            *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
            this->GetParamData(), this->GetDomainData());
        auto& edc = GetIntegratorDataContainer().GetElementDataContainer();

        dwrc.GetWeightIDC().InitializeEDC(pde.GetUpdateFlags(),
            dwrc.GetWeightSTH(), element_weight, this->GetParamData(),
            dwrc.GetWeightData());
        auto& edc_weight = dwrc.GetElementWeight();

        // we want to integrate the face-terms only once, so
        // we store the values on each face in this map
        // and distribute it at the end to the adjacent elements.
        typename std::map<typename dealii::Triangulation<dim>::face_iterator,std::vector<double> >
        face_integrals;
        // initialize the map with a big value to make sure
        // that we take notice if we forget to add a face
        // during the error estimation process
        auto element_it = element[0];
        std::vector<double> face_init(n_error_comps,-1e20);
        for (; element_it != endc[0]; element_it++)
        {
          for (unsigned int face_no=0;face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
          {
            face_integrals[element_it->face(face_no)] = face_init;
          }
        }

//        bool need_faces = pde.HasFaces();
//        bool need_interfaces = pde.HasInterfaces();
//        std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
//        bool need_boundary_integrals = (boundary_equation_colors.size() > 0);

        GetIntegratorDataContainer().InitializeFDC(dwrc.GetWeightIDC().GetFaceQuad(),
            pde.GetFaceUpdateFlags(),
            *(pde.GetBaseProblem().GetSpaceTimeHandler()),
            element,
            this->GetParamData(),
            this->GetDomainData(),
            true);
        auto & fdc = GetIntegratorDataContainer().GetFaceDataContainer();

        dwrc.GetWeightIDC().InitializeFDC(pde.GetFaceUpdateFlags(),
            dwrc.GetWeightSTH(),
            element_weight,
            this->GetParamData(),
            dwrc.GetWeightData(),
            true);

        for (unsigned int element_index = 0; element[0] != endc[0];
            element[0]++, element_index++)
        {
          for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
          {
            if (element[dh] == endc[dh])
            {
              throw DOpEException(
                  "Elementnumbers in DoFHandlers are not matching!",
                  "Integrator::ComputeRefinementIndicators");
            }
          }
          for (unsigned int dh = 0; dh < dof_handler_weight.size(); dh++)
          {
            if (element_weight[dh] == endc_high[dh])
            {
              throw DOpEException(
                  "Elementnumbers in DoFHandlers are not matching!",
                  "Integrator::ComputeRefinementIndicators");
            }
          }
          element_sum.clear();
          element_sum.resize(n_error_comps, 0);

          edc.ReInit();
          edc_weight.ReInit();

          //first the element-residual
          pde.ElementErrorContribution(edc, dwrc, element_sum, 1.);
          for(unsigned int l =0; l < n_error_comps; l ++)
          {
            dwrc.GetErrorIndicators(l)(element_index) = element_sum[l];
          }
          //dwrc.GetPrimalErrorIndicators()(element_index) = element_sum[0];
          //dwrc.GetDualErrorIndicators()(element_index) = element_sum[1];
          element_sum.clear();
          element_sum.resize(n_error_comps, 0);


          //Now to the face terms. We compute them only once for each face and distribute the
          //afterwards. We choose always to work from the coarser element, if both neigbors of the
          //face are on the same level, we pick the one with the lower index
          for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
          {
            auto face_it = element[0]->face(face);

            //check if the face lies at a boundary
            if(face_it->at_boundary())
            {
              fdc.ReInit(face);
              dwrc.GetFaceWeight().ReInit(face);
              pde.BoundaryErrorContribution(fdc, dwrc, element_sum, 1.);

              Assert (face_integrals.find (element[0]->face(face)) != face_integrals.end(),
                  ExcInternalError());
              Assert (face_integrals[element[0]->face(face)] == face_init,
                  ExcInternalError());

              face_integrals[element[0]->face(face)] = element_sum;
              element_sum.clear();
              element_sum.resize(n_error_comps,0.);
            }
            else
            {
              //There exist now 3 different scenarios, given the actual element and face:
              // The neighbour behind this face is [ more | as much | less] refined
              // than/as the actual element. We have to distinguish here only between the case 1
              // and the other two, because these will be distinguished in in the FaceDataContainer.
              if (element[0]->neighbor(face)->has_children())
              {
                //first: neighbour is finer
                std::vector<double> sum(n_error_comps,0.);
                for (unsigned int subface_no=0;
                    subface_no < element[0]->face(face)->n_children();
                    ++subface_no)
                {
                  //TODO Now here we have to initialise the subface_values on the
                  // actual element and then the facevalues of the neighbours
                  fdc.ReInit(face, subface_no);
                  fdc.ReInitNbr();
                  dwrc.GetFaceWeight().ReInit(face, subface_no);

                  pde.FaceErrorContribution(fdc, dwrc, element_sum, 1.);
                  for(unsigned int l =0; l < n_error_comps; l ++)
                  {
                    sum[l] += element_sum[l];
                  }
                  face_integrals[element[0]->neighbor_child_on_subface(face, subface_no)
                  ->face(element[0]->neighbor_of_neighbor(face))] = element_sum;
                  element_sum.clear();
                  element_sum.resize(n_error_comps,0);
                }

                Assert (face_integrals.find (element[0]->face(face)) != face_integrals.end(),
                    ExcInternalError());
                Assert (face_integrals[element[0]->face(face)] == face_init,
                    ExcInternalError());

                face_integrals[element[0]->face(face)] = sum;
                element_sum.clear();
                element_sum.resize(n_error_comps,0);
              }
              else
              {
                // either neighbor is as fine as this element or
                // it is coarser
                Assert(element[0]->neighbor(face)->level() <= element[0]->level(),ExcInternalError());
                //now we work always from the coarser element. if both elements
                //are on the same level, we pick the one with the lower index
                if(element[0]->level() == element[0]->neighbor(face)->level()
                    && element[0]->index() < element[0]->neighbor(face)->index())
                {
                  fdc.ReInit(face);
                  fdc.ReInitNbr();
                  dwrc.GetFaceWeight().ReInit(face);

                  pde.FaceErrorContribution(fdc, dwrc, element_sum, 1.);
                  Assert (face_integrals.find (element[0]->face(face)) != face_integrals.end(),
                      ExcInternalError());
                  Assert (face_integrals[element[0]->face(face)] == face_init,
                      ExcInternalError());

                  face_integrals[element[0]->face(face)] = element_sum;
                  element_sum.clear();
                  element_sum.resize(n_error_comps,0);
                }
              }
            }
          }                  //endfor faces
          for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
          {
            element[dh]++;
          }
          for (unsigned int dh = 0; dh < dof_handler_weight.size(); dh++)
          {
            element_weight[dh]++;
          }
        }                  //endfor element
        //now we have to incorporate the face and boundary_values
        //into
        unsigned int present_element = 0;
        element =
        pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
        for (;
            element[0] !=endc[0]; element[0]++, ++present_element)
        {
          for (unsigned int face_no = 0;
              face_no < GeometryInfo<dim>::faces_per_cell; ++face_no)
          {
            Assert(
                face_integrals.find(element[0]->face(face_no)) != face_integrals.end(),
                ExcInternalError());
            if(element[0]->face(face_no)->at_boundary())
            {
              for(unsigned int l =0; l < n_error_comps; l ++)
              {
                dwrc.GetErrorIndicators(l)(present_element) +=
                  face_integrals[element[0]->face(face_no)][l];
              }
//              dwrc.GetPrimalErrorIndicators()(present_element) +=
//              face_integrals[element[0]->face(face_no)][0];
//              dwrc.GetDualErrorIndicators()(present_element) +=
//              face_integrals[element[0]->face(face_no)][1];
            }
            else
            {
              for(unsigned int l =0; l < n_error_comps; l ++)
              {
                dwrc.GetErrorIndicators(l)(present_element) +=
                  0.5*face_integrals[element[0]->face(face_no)][l];
              }
//              dwrc.GetPrimalErrorIndicators()(present_element) +=
//              0.5 * face_integrals[element[0]->face(face_no)][0];
//              dwrc.GetDualErrorIndicators()(present_element) +=
//              0.5 * face_integrals[element[0]->face(face_no)][1];
            }

          }
        }
      }
  /*******************************************************************************************/

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<typename PROBLEM>
      void
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeRefinementIndicators(
          PROBLEM& pde,
          ResidualErrorContainer<VECTOR>& dwrc)
      {
        //for primal and dual part of the error
        std::vector<double> element_sum(2, 0);
        element_sum.resize(2, 0);
        // Begin integration
        const auto& dof_handler =
            pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
        auto element =
            pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
        auto endc =
            pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();

        {
          //Add Weights
          auto wd = dwrc.GetWeightData().begin();
          auto wend = dwrc.GetWeightData().end();
          for (; wd != wend; wd++)
          {
            AddDomainData(wd->first, wd->second);
          }
        }

        // Generate the data containers. Notice that we use the quadrature
        // formula from the higher order idc!.
        GetIntegratorDataContainer().InitializeEDC(pde.GetUpdateFlags(),
            *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
            this->GetParamData(), this->GetDomainData());
        auto& edc = GetIntegratorDataContainer().GetElementDataContainer();

        //we want to integrate the face-terms only once
        typename std::map<typename dealii::Triangulation<dim>::face_iterator,std::vector<double> >
        face_integrals;
        //initialize the map
        auto element_it = element[0];
        std::vector<double> face_init(2,-1e20);
        for (; element_it != endc[0]; element_it++)
        {
          for (unsigned int face_no=0;face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
          {
            face_integrals[element_it->face(face_no)] = face_init;
          }
        }

//        bool need_faces = pde.HasFaces();
//        bool need_interfaces = pde.HasInterfaces();
//        std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
//        bool need_boundary_integrals = (boundary_equation_colors.size() > 0);

        GetIntegratorDataContainer().InitializeFDC(
            pde.GetFaceUpdateFlags(),
            *(pde.GetBaseProblem().GetSpaceTimeHandler()),
            element,
            this->GetParamData(),
            this->GetDomainData(),
            true);
        auto & fdc = GetIntegratorDataContainer().GetFaceDataContainer();

        for (unsigned int element_index = 0; element[0] != endc[0];
            element[0]++, element_index++)
        {
          for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
          {
            if (element[dh] == endc[dh])
            {
              throw DOpEException(
                  "Elementnumbers in DoFHandlers are not matching!",
                  "Integrator::ComputeRefinementIndicators");
            }
          }

          element_sum.clear();
          element_sum.resize(2, 0);

          edc.ReInit();
          dwrc.InitElement(element[0]->diameter());
          //first the element-residual
          pde.ElementErrorContribution(edc, dwrc, element_sum, 1.);
          dwrc.GetPrimalErrorIndicators()(element_index) = element_sum[0];
          dwrc.GetDualErrorIndicators()(element_index) = element_sum[1];
          element_sum.clear();
          element_sum.resize(2, 0);
          //Now to the face terms. We compute them only once for each face and distribute the
          //afterwards. We choose always to work from the coarser element, if both neigbors of the
          //face are on the same level, we pick the one with the lower index

          for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
          {
            auto face_it = element[0]->face(face);

            //check if the face lies at a boundary
            if(face_it->at_boundary())
            {
              fdc.ReInit(face);
#if deal_II_dimension > 1
              dwrc.InitFace(element[0]->face(face)->diameter());
#else 
              dwrc.InitFace(element[0]->diameter());
#endif   
              pde.BoundaryErrorContribution(fdc, dwrc, element_sum, 1.);

              Assert (face_integrals.find (element[0]->face(face)) != face_integrals.end(),
                  ExcInternalError());
              Assert (face_integrals[element[0]->face(face)] == face_init,
                  ExcInternalError());
              face_integrals[element[0]->face(face)] = element_sum;
              element_sum.clear();
              element_sum.resize(2,0.);
            }
            else
            {
              //There exist now 3 different scenarios, given the actual element and face:
              // The neighbour behind this face is [ more | as much | less] refined
              // than/as the actual element. We have to distinguish here only between the case 1
              // and the other two, because these will be distinguished in in the FaceDataContainer.
              if (element[0]->neighbor(face)->has_children())
              {
                //first: neighbour is finer
                std::vector<double> sum(2,0.);
                for (unsigned int subface_no=0;
                    subface_no < element[0]->face(face)->n_children();
                    ++subface_no)
                {
                  //TODO Now here we have to initialise the subface_values on the
                  // actual element and then the facevalues of the neighbours
                  fdc.ReInit(face, subface_no);
                  fdc.ReInitNbr();
#if deal_II_dimension > 1
                  dwrc.InitFace(element[0]->face(face)->diameter());
#else 
                  dwrc.InitFace(element[0]->diameter());
#endif   

                  pde.FaceErrorContribution(fdc, dwrc, element_sum, 1.);
                  sum[0]= element_sum[0];
                  sum[1]= element_sum[1];
                  element_sum.clear();
                  element_sum.resize(2,0);
                  face_integrals[element[0]->neighbor_child_on_subface(face, subface_no)
                  ->face(element[0]->neighbor_of_neighbor(face))] = element_sum;
                  element_sum.clear();
                  element_sum.resize(2,0.);
                }

                Assert (face_integrals.find (element[0]->face(face)) != face_integrals.end(),
                    ExcInternalError());
                Assert (face_integrals[element[0]->face(face)] == face_init,
                    ExcInternalError());

                face_integrals[element[0]->face(face)] = sum;

              }
              else
              {
                // either neighbor is as fine as this element or
                // it is coarser
                Assert(element[0]->neighbor(face)->level() <= element[0]->level(),ExcInternalError());
                //now we work always from the coarser element. if both elements
                //are on the same level, we pick the one with the lower index
                if(element[0]->level() == element[0]->neighbor(face)->level()
                    && element[0]->index() < element[0]->neighbor(face)->index())
                {
                  fdc.ReInit(face);
                  fdc.ReInitNbr();
#if deal_II_dimension > 1
                  dwrc.InitFace(element[0]->face(face)->diameter());
#else 
                  dwrc.InitFace(element[0]->diameter());
#endif   
         
                  pde.FaceErrorContribution(fdc, dwrc, element_sum, 1.);
                  Assert (face_integrals.find (element[0]->face(face)) != face_integrals.end(),
                      ExcInternalError());
                  Assert (face_integrals[element[0]->face(face)] == face_init,
                      ExcInternalError());

                  face_integrals[element[0]->face(face)] = element_sum;
                  element_sum.clear();
                  element_sum.resize(2,0);
                }
              }
            }
          }                  //endfor faces
//          }//end else

          for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
          {
            element[dh]++;
          }
        }                  //endfor element
        //now we have to incorporate the face and boundary_values
        //into
        unsigned int present_element = 0;
        element =
        pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
        for (;
            element[0] !=endc[0]; element[0]++, ++present_element)
        for (unsigned int face_no = 0;
            face_no < GeometryInfo<dim>::faces_per_cell; ++face_no)
        {
          Assert(
              face_integrals.find(element[0]->face(face_no)) != face_integrals.end(),
              ExcInternalError());
          dwrc.GetPrimalErrorIndicators()(present_element) +=
          0.5 * face_integrals[element[0]->face(face_no)][0];
          dwrc.GetDualErrorIndicators()(present_element) +=
          0.5 * face_integrals[element[0]->face(face_no)][1];
        }
        {
          //Remove Weights
          auto wd = dwrc.GetWeightData().begin();
          auto wend = dwrc.GetWeightData().end();
          for (; wd != wend; wd++)
          {
            DeleteDomainData(wd->first);
          }
        }
      }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    INTEGRATORDATACONT&
    Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::GetIntegratorDataContainer() const
    {
      return _idc1;
    }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    INTEGRATORDATACONT&
    Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::GetIntegratorDataContainerFunc() const
    {
      return _idc2;
    }

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

  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
      int dim>
    template<template<int, int> class DH>
      void
      Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::InterpolateBoundaryValues(
          const DOpEWrapper::Mapping<dim, DH>& mapping,
          const DOpEWrapper::DoFHandler<dim, DH>* dof_handler,
          const unsigned int color, const dealii::Function<dim>& function,
          std::map<unsigned int, SCALAR>& boundary_values,
          const std::vector<bool>& comp_mask) const
      {
    //TODO: mapping[0] is a workaround, as deal does not support interpolate
    // boundary_values with a mapping collection at this point.
        dealii::VectorTools::interpolate_boundary_values(mapping[0],
            dof_handler->GetDEALDoFHandler(), color, function,
            boundary_values, comp_mask);
      }


}
#endif