timedofhandler.h

Go to the documentation of this file.

#ifndef TIMEDOFHANDLER_H_
#define TIMEDOFHANDLER_H_

//DOpE
#include <timeiterator.h>
#include <dopeexception.h>
#include <versionscheck.h>

//deal.ii
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_update_flags.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/base/point.h>

//c++
#include <vector>


using namespace dealii;

namespace DOpE
{

  class TimeDoFHandler : public dealii::DoFHandler<1>
  {
    public:
      //Constructors
      TimeDoFHandler()
      {
        fe_ = NULL;
        times_.resize(1, 0.);
        dofs_per_element_ = 1;
        need_delete_ = false;
        initialized_ = false;
      }

      TimeDoFHandler(const Triangulation<1> & tria,
          const dealii::FiniteElement<1>& fe)
          : dealii::DoFHandler<1>(tria)
      {
        fe_ = &fe;
        this->distribute_dofs();
        need_delete_ = false;
      }

      TimeDoFHandler(const Triangulation<1>& tria)
          : dealii::DoFHandler<1>(tria)
      {
        fe_ = new dealii::FE_Q<1>(1);
        this->distribute_dofs();
        need_delete_ = true;
      }

      //Destructor
      ~TimeDoFHandler()
      {
        this->clear();

        if (need_delete_)
        {
          delete fe_;
        }
        fe_ = NULL;
      }

      void
      distribute_dofs()
      {
        if (fe_ != NULL)
        {
          dealii::DoFHandler<1>::distribute_dofs(*fe_);
          //make sure that the dofs are numbered 'downstream' (referring to the time variable!)

#if DEAL_II_VERSION_GTE(7,3,0)
           dealii::DoFRenumbering::downstream<dealii::DoFHandler<1> >(*this,
               dealii::Point<1>(1.), true);
  #else
           dealii::DoFRenumbering::downstream<dealii::DoFHandler<1>, 1>(*this,
               dealii::Point<1>(1.), true);
  #endif
          find_ends();
          compute_times();
          initialized_ = true;
          dofs_per_element_ = this->get_fe().dofs_per_cell;
        }
      }

      unsigned int
      GetNbrOfIntervals() const
      {
        if (initialized_)
          return this->get_tria().n_active_cells();
        else
          return 0;
      }

      unsigned int
      GetNbrOfDoFs() const
      {
        if (initialized_)
          return this->n_dofs();
        else
          return 1;
      }

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

      void
      GetTimes(const TimeIterator & interval,
          std::vector<double>& local_times) const
      {
        assert(local_times.size() == GetLocalNbrOfDoFs());

        std::vector<unsigned int> global_dof_indices(GetLocalNbrOfDoFs());
        interval.get_time_dof_indices(global_dof_indices);
        for (unsigned int i = 0; i < GetLocalNbrOfDoFs(); ++i)
        {
          local_times[i] = times_[global_dof_indices[i]];
        }
      }

      double
      GetTime(unsigned int timestep)
      {
        assert(timestep < times_.size());
        return times_[timestep];
      }

      TimeIterator
      first_interval() const
      {
        assert(initialized_);
        return first_interval_;
      }

      TimeIterator
      before_first_interval() const
      {
        assert(initialized_);
        return before_first_interval_;
      }

      TimeIterator
      last_interval() const
      {
        assert(initialized_);
        return last_interval_;
      }

      TimeIterator
      after_last_interval() const
      {
        assert(initialized_);
        return after_last_interval_;
      }

      unsigned int
      GetLocalNbrOfDoFs() const
      {
        return dofs_per_element_;
      }

    private:
      using dealii::DoFHandler<1>::distribute_dofs;   //Remove the clang warning that the function below overloads the distribute_dofs
      void
      find_ends()
      {
        DoFHandler<1>::active_cell_iterator element = this->begin_active();
#if DEAL_II_VERSION_GTE(8,3,0)
        while (element->face(0)->boundary_id() != 0)
#else
        while (element->face(0)->boundary_indicator() != 0)
#endif
        {
          element = element->neighbor(0);
        }
        first_interval_.Initialize(element, 0);
        before_first_interval_.Initialize(element, -2);

        element = this->begin_active();
#if DEAL_II_VERSION_GTE(8,3,0)
        while (element->face(1)->boundary_id() != 1)
#else
        while (element->face(1)->boundary_indicator() != 1)
#endif
        {
          element = element->neighbor(1);
        }
        assert(static_cast<int>(this->get_tria().n_active_cells() - 1) >= 0);
        last_interval_.Initialize(element,
            static_cast<int>(this->get_tria().n_active_cells() - 1));after_last_interval_
        .Initialize(element, -1);
      }

      void
      compute_times()
      {
        assert(this->get_fe().has_support_points());
        //get the right length
        times_.resize(this->n_dofs());

        //geht the support points and build a quadrature rule from it
        dealii::Quadrature<1> quadrature_formula(
            this->get_fe().get_unit_support_points());

        //after that, build a fevalues object. We need only the quadrature points!
        //typename dealii::FEValues<1> fe_values(this->get_fe(),
        dealii::FEValues<1> fe_values(this->get_fe(), quadrature_formula,
            dealii::UpdateFlags(update_q_points));

        const unsigned int n_q_points = quadrature_formula.size();
        std::vector<unsigned int> global_dof_indices(
            this->get_fe().dofs_per_cell);

        //the usual loops. Go over all elements, in every element go through the quad
        //points and store thei position in times_
        //typename DoFHandler<1>::active_cell_iterator element =
        DoFHandler<1>::active_cell_iterator element = this->begin_active(), endc =
            this->end();
        for (; element != endc; ++element)
        {
          fe_values.reinit(element);
          element->get_dof_indices(global_dof_indices);
          for (unsigned int i = 0; i < n_q_points; ++i)
          {
            times_[global_dof_indices[i]] =
                fe_values.get_quadrature_points()[i](0);
          }
        }
      }
      //member variables

      TimeIterator first_interval_, last_interval_, before_first_interval_,
          after_last_interval_;

      //FIXME
      //Dummy FE which is used, if no fe is specified from the user.
      //At the moment, this is necessary, because we use an DoFCellAccessor
      //-iterator to go through our timegrid. Is this ok?
      const dealii::FiniteElement<1>* fe_;

      std::vector<double> times_;
      unsigned int dofs_per_element_;
      bool need_delete_, initialized_;
  };

}      //end of namespace

#endif /* TIMEDOFHANDLER_H_ */