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)
           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:
      void
      find_ends()
      {
        DoFHandler<1>::active_cell_iterator element = this->begin_active();
        while (element->face(0)->boundary_indicator() != 0)
        {
          element = element->neighbor(0);
        }
        _first_interval.Initialize(element, 0);
        _before_first_interval.Initialize(element, -2);

        element = this->begin_active();
        while (element->face(1)->boundary_indicator() != 1)
        {
          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_ */