dof_info.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2015 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------


#ifndef dealii__matrix_free_dof_info_h
#define dealii__matrix_free_dof_info_h


#include <deal.II/base/exceptions.h>
#include <deal.II/base/vectorization.h>
#include <deal.II/base/partitioner.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/matrix_free/helper_functions.h>

#include <deal.II/base/std_cxx11/array.h>

#include <memory>


DEAL_II_NAMESPACE_OPEN

namespace internal
{
  namespace MatrixFreeFunctions
  {
    struct DoFInfo
    {
      DoFInfo ();

      DoFInfo (const DoFInfo &dof_info);

      void clear ();


      const unsigned int *begin_indices (const unsigned int row) const;

      const unsigned int *end_indices (const unsigned int row) const;

      unsigned int row_length_indices (const unsigned int row) const;

      const std::pair<unsigned short,unsigned short> *
      begin_indicators (const unsigned int row) const;

      const std::pair<unsigned short,unsigned short> *
      end_indicators (const unsigned int row) const;

      unsigned int row_length_indicators (const unsigned int row) const;

      const unsigned int *begin_indices_plain (const unsigned int row) const;

      const unsigned int *end_indices_plain (const unsigned int row) const;

      unsigned int fe_index_from_degree (const unsigned int fe_degree) const;


      unsigned int
      fe_index_from_dofs_per_cell (const unsigned int dofs_per_cell) const;

      void read_dof_indices (const std::vector<types::global_dof_index> &local_indices,
                             const std::vector<unsigned int> &lexicographic_inv,
                             const ConstraintMatrix          &constraints,
                             const unsigned int               cell_number,
                             ConstraintValues<double> &constraint_values,
                             bool                            &cell_at_boundary);

      void assign_ghosts(const std::vector<unsigned int> &boundary_cells);

      void compute_renumber_serial (const std::vector<unsigned int> &boundary_cells,
                                    const SizeInfo                  &size_info,
                                    std::vector<unsigned int>       &renumbering);

      void compute_renumber_hp_serial (SizeInfo                  &size_info,
                                       std::vector<unsigned int> &renumbering,
                                       std::vector<unsigned int> &irregular_cells);

      void compute_renumber_parallel (const std::vector<unsigned int> &boundary_cells,
                                      SizeInfo                        &size_info,
                                      std::vector<unsigned int>       &renumbering);

      void reorder_cells (const SizeInfo                  &size_info,
                          const std::vector<unsigned int> &renumbering,
                          const std::vector<unsigned int> &constraint_pool_row_index,
                          const std::vector<unsigned int> &irregular_cells,
                          const unsigned int               vectorization_length);

      void guess_block_size (const SizeInfo &size_info,
                             TaskInfo       &task_info);

      void
      make_thread_graph_partition_color (SizeInfo                  &size_info,
                                         TaskInfo                  &task_info,
                                         std::vector<unsigned int> &renumbering,
                                         std::vector<unsigned int> &irregular_cells,
                                         const bool                 hp_bool);

      void
      make_thread_graph_partition_partition (SizeInfo                  &size_info,
                                             TaskInfo                  &task_info,
                                             std::vector<unsigned int> &renumbering,
                                             std::vector<unsigned int> &irregular_cells,
                                             const bool                 hp_bool);

      void
      make_connectivity_graph (const SizeInfo                  &size_info,
                               const TaskInfo                  &task_info,
                               const std::vector<unsigned int> &renumbering,
                               const std::vector<unsigned int> &irregular_cells,
                               const bool                       do_blocking,
                               DynamicSparsityPattern &connectivity) const;

      void renumber_dofs (std::vector<types::global_dof_index> &renumbering);

      std::size_t memory_consumption() const;

      template <typename StreamType>
      void print_memory_consumption(StreamType     &out,
                                    const SizeInfo &size_info) const;

      template <typename Number>
      void print (const std::vector<Number>       &constraint_pool_data,
                  const std::vector<unsigned int> &constraint_pool_row_index,
                  std::ostream                    &out) const;

      std::vector<std_cxx11::array<unsigned int, 3> > row_starts;

      std::vector<unsigned int> dof_indices;

      std::vector<std::pair<unsigned short,unsigned short> > constraint_indicator;

      std_cxx11::shared_ptr<const Utilities::MPI::Partitioner> vector_partitioner;

      std::vector<unsigned int> constrained_dofs;

      std::vector<unsigned int> row_starts_plain_indices;

      std::vector<unsigned int> plain_dof_indices;

      unsigned int dimension;

      unsigned int n_components;

      std::vector<unsigned int> dofs_per_cell;

      std::vector<unsigned int> dofs_per_face;

      bool store_plain_indices;

      std::vector<unsigned int> cell_active_fe_index;

      unsigned int max_fe_index;

      std::vector<std::pair<unsigned int,unsigned int> > fe_index_conversion;

      std::vector<types::global_dof_index> ghost_dofs;
    };


    /*----------------------- Inline functions ----------------------------------*/

#ifndef DOXYGEN

    inline
    const unsigned int *
    DoFInfo::begin_indices (const unsigned int row) const
    {
      AssertIndexRange (row, row_starts.size()-1);
      const unsigned int index = row_starts[row][0];
      AssertIndexRange(index, dof_indices.size()+1);
      return dof_indices.empty() ? 0 : &dof_indices[0] + index;
    }



    inline
    const unsigned int *
    DoFInfo::end_indices (const unsigned int row) const
    {
      AssertIndexRange (row, row_starts.size()-1);
      const unsigned int index = row_starts[row+1][0];
      AssertIndexRange(index, dof_indices.size()+1);
      return dof_indices.empty() ? 0 : &dof_indices[0] + index;
    }



    inline
    unsigned int
    DoFInfo::row_length_indices (const unsigned int row) const
    {
      AssertIndexRange (row, row_starts.size()-1);
      return (row_starts[row+1][0] - row_starts[row][0]);
    }



    inline
    const std::pair<unsigned short,unsigned short> *
    DoFInfo::begin_indicators (const unsigned int row) const
    {
      AssertIndexRange (row, row_starts.size()-1);
      const unsigned int index = row_starts[row][1];
      AssertIndexRange (index, constraint_indicator.size()+1);
      return constraint_indicator.empty() ? 0 : &constraint_indicator[0] + index;
    }



    inline
    const std::pair<unsigned short,unsigned short> *
    DoFInfo::end_indicators (const unsigned int row) const
    {
      AssertIndexRange (row, row_starts.size()-1);
      const unsigned int index = row_starts[row+1][1];
      AssertIndexRange (index, constraint_indicator.size()+1);
      return constraint_indicator.empty() ? 0 : &constraint_indicator[0] + index;
    }



    inline
    unsigned int
    DoFInfo::row_length_indicators (const unsigned int row) const
    {
      AssertIndexRange (row, row_starts.size()-1);
      return (row_starts[row+1][1] - row_starts[row][1]);
    }



    inline
    const unsigned int *
    DoFInfo::begin_indices_plain (const unsigned int row) const
    {
      // if we have no constraints, should take the data from dof_indices
      if (row_length_indicators(row) == 0)
        {
          Assert (row_starts_plain_indices[row]==numbers::invalid_unsigned_int,
                  ExcInternalError());
          return begin_indices(row);
        }
      else
        {
          AssertDimension (row_starts.size(), row_starts_plain_indices.size());
          const unsigned int index = row_starts_plain_indices[row];
          AssertIndexRange(index, plain_dof_indices.size()+1);
          return plain_dof_indices.empty() ? 0 : &plain_dof_indices[0] + index;
        }
    }



    inline
    const unsigned int *
    DoFInfo::end_indices_plain (const unsigned int row) const
    {
      return begin_indices_plain(row) +
             dofs_per_cell[(cell_active_fe_index.size()==0)?
                           0:cell_active_fe_index[row]];
    }



    inline
    unsigned int
    DoFInfo::fe_index_from_degree (const unsigned int fe_degree) const
    {
      const unsigned int n_indices = fe_index_conversion.size();
      for (unsigned int i=0; i<n_indices; ++i)
        if (fe_index_conversion[i].first == fe_degree)
          return i;
      return n_indices;
    }



    inline
    unsigned int
    DoFInfo::fe_index_from_dofs_per_cell (const unsigned int dofs_per_cell) const
    {
      for (unsigned int i=0; i<fe_index_conversion.size(); ++i)
        if (fe_index_conversion[i].second == dofs_per_cell)
          return i;
      return 0;
    }

  } // end of namespace MatrixFreeFunctions
} // end of namespace internal

#endif  // ifndef DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif