<commit_before>// ---------------------------------------------------------------------
//
// Copyright (C) 2003 - 2018 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.
//
// ---------------------------------------------------------------------

#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/dof_level.h>
#include <deal.II/hp/dof_faces.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_handler_policy.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_levels.h>
#include <deal.II/grid/tria.h>
#include <deal.II/fe/fe.h>
#include <deal.II/distributed/shared_tria.h>
#include <deal.II/distributed/tria.h>

#include <boost/serialization/array.hpp>

#include <set>
#include <algorithm>
#include <functional>

DEAL_II_NAMESPACE_OPEN

// The following is necessary for compilation under Visual Studio which is unable to correctly
// distinguish between dealii::DoFHandler and dealii::hp::DoFHandler.
// Plus it makes code in dof_handler.cc easier to read.
// Requires C++11 support which is in Visual Studio 2013 and newer.
#if defined(_MSC_VER) && (_MSC_VER >= 1800)
template <int dim, int spacedim> using HpDoFHandler = ::dealii::hp::DoFHandler<dim, spacedim>;
#else
// When using older Visual Studio or a different compiler just fall back.
#define HpDoFHandler DoFHandler
#endif

namespace parallel
{
  namespace distributed
  {
    template <int, int> class Triangulation;
  }
}




namespace internal
{
  namespace hp
  {
    namespace DoFHandler
    {
      // access class dealii::hp::DoFHandler instead of namespace
      // internal::hp::DoFHandler, etc
      using dealii::hp::DoFHandler;

      /**
       * A class with the same purpose as the similarly named class of the
       * Triangulation class. See there for more information.
       */
      struct Implementation
      {
        /**
         * Do that part of reserving space that pertains to releasing
         * the previously used memory.
         */
        template <int dim, int spacedim>
        static
        void
        reserve_space_release_space (DoFHandler<dim,spacedim> &dof_handler)
        {
          // Release all space except the active_fe_indices field
          // which we have to back up before
          {
            std::vector<std::vector<DoFLevel::active_fe_index_type> >
            active_fe_backup(dof_handler.levels.size ());
            for (unsigned int level = 0; level<dof_handler.levels.size (); ++level)
              active_fe_backup[level] = std::move(dof_handler.levels[level]->active_fe_indices);

            // delete all levels and set them up newly, since vectors
            // are troublesome if you want to change their size
            dof_handler.clear_space ();

            for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
              {
                dof_handler.levels.emplace_back (new internal::hp::DoFLevel);
                dof_handler.levels[level]->active_fe_indices = std::move(active_fe_backup[level]);
              }

            if (dim > 1)
              dof_handler.faces.reset (new internal::hp::DoFIndicesOnFaces<dim>);
          }
        }



        /**
         * Do that part of reserving space that pertains to vertices,
         * since this is the same in all space dimensions.
         */
        template <int dim, int spacedim>
        static
        void
        reserve_space_vertices (DoFHandler<dim,spacedim> &dof_handler)
        {
          // The final step in all of the reserve_space() functions is to set
          // up vertex dof information. since vertices are sequentially
          // numbered, what we do first is to set up an array in which
          // we record whether a vertex is associated with any of the
          // given fe's, by setting a bit. in a later step, we then
          // actually allocate memory for the required dofs
          //
          // in the following, we only need to consider vertices that are
          // adjacent to either a locally owned or a ghost cell; we never
          // store anything on vertices that are only surrounded by
          // artificial cells. so figure out that subset of vertices
          // first
          std::vector<bool> locally_used_vertices (dof_handler.tria->n_vertices(),
                                                   false);
          for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
               cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
            if (! cell->is_artificial())
              for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
                locally_used_vertices[cell->vertex_index(v)] = true;

          std::vector<std::vector<bool> >
          vertex_fe_association (dof_handler.finite_elements->size(),
                                 std::vector<bool> (dof_handler.tria->n_vertices(), false));

          for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
               cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
            if (! cell->is_artificial())
              for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
                vertex_fe_association[cell->active_fe_index()][cell->vertex_index(v)]
                  = true;

          // in debug mode, make sure that each vertex is associated
          // with at least one fe (note that except for unused
          // vertices, all vertices are actually active). this is of
          // course only true for vertices that are part of either
          // ghost or locally owned cells
#ifdef DEBUG
          for (unsigned int v=0; v<dof_handler.tria->n_vertices(); ++v)
            if (locally_used_vertices[v] == true)
              if (dof_handler.tria->vertex_used(v) == true)
                {
                  unsigned int fe=0;
                  for (; fe<dof_handler.finite_elements->size(); ++fe)
                    if (vertex_fe_association[fe][v] == true)
                      break;
                  Assert (fe != dof_handler.finite_elements->size(), ExcInternalError());
                }
#endif

          // next count how much memory we actually need. for each
          // vertex, we need one slot per fe to store the fe_index,
          // plus dofs_per_vertex for this fe. in addition, we need
          // one slot as the end marker for the fe_indices. at the
          // same time already fill the vertex_dof_offsets field
          dof_handler.vertex_dof_offsets.resize (dof_handler.tria->n_vertices(),
                                                 numbers::invalid_unsigned_int);

          unsigned int vertex_slots_needed = 0;
          for (unsigned int v=0; v<dof_handler.tria->n_vertices(); ++v)
            if (dof_handler.tria->vertex_used(v) == true)
              if (locally_used_vertices[v] == true)
                {
                  dof_handler.vertex_dof_offsets[v] = vertex_slots_needed;

                  for (unsigned int fe=0; fe<dof_handler.finite_elements->size(); ++fe)
                    if (vertex_fe_association[fe][v] == true)
                      vertex_slots_needed += dof_handler.get_fe(fe).dofs_per_vertex + 1;

                  // don't forget the end_marker:
                  ++vertex_slots_needed;
                }

          // now allocate the space we have determined we need, and
          // set up the linked lists for each of the vertices
          dof_handler.vertex_dofs.resize (vertex_slots_needed,
                                          numbers::invalid_dof_index);
          for (unsigned int v=0; v<dof_handler.tria->n_vertices(); ++v)
            if (dof_handler.tria->vertex_used(v) == true)
              if (locally_used_vertices[v] == true)
                {
                  unsigned int current_index = dof_handler.vertex_dof_offsets[v];
                  for (unsigned int fe=0; fe<dof_handler.finite_elements->size(); ++fe)
                    if (vertex_fe_association[fe][v] == true)
                      {
                        // if this vertex uses this fe, then set the
                        // fe_index and move the pointer ahead
                        dof_handler.vertex_dofs[current_index] = fe;
                        current_index += dof_handler.get_fe(fe).dofs_per_vertex + 1;
                      }
                  // finally place the end marker
                  dof_handler.vertex_dofs[current_index] = numbers::invalid_dof_index;
                }
        }


        /**
         * Do that part of reserving space that pertains to cells,
         * since this is the same in all space dimensions.
         */
        template <int dim, int spacedim>
        static
        void
        reserve_space_cells (DoFHandler<dim,spacedim> &dof_handler)
        {
          // count how much space we need on each level for the cell
          // dofs and set the dof_*_offsets data. initially set the
          // latter to an invalid index, and only later set it to
          // something reasonable for active dof_handler.cells
          //
          // note that for dof_handler.cells, the situation is simpler
          // than for other (lower dimensional) objects since exactly
          // one finite element is used for it
          for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
            {
              dof_handler.levels[level]->dof_offsets
                = std::vector<DoFLevel::offset_type> (dof_handler.tria->n_raw_cells(level),
                                                      (DoFLevel::offset_type)(-1));
              dof_handler.levels[level]->cell_cache_offsets
                = std::vector<DoFLevel::offset_type> (dof_handler.tria->n_raw_cells(level),
                                                      (DoFLevel::offset_type)(-1));

              types::global_dof_index next_free_dof = 0;
              types::global_dof_index cache_size = 0;
              typename HpDoFHandler<dim,spacedim>::active_cell_iterator
              cell=dof_handler.begin_active(level),
              endc=dof_handler.end_active(level);
              for (; cell!=endc; ++cell)
                if (!cell->has_children()
                    &&
                    !cell->is_artificial())
                  {
                    dof_handler.levels[level]->dof_offsets[cell->index()] = next_free_dof;
                    next_free_dof += cell->get_fe().template n_dofs_per_object<dim>();

                    dof_handler.levels[level]->cell_cache_offsets[cell->index()] = cache_size;
                    cache_size += cell->get_fe().dofs_per_cell;
                  }

              dof_handler.levels[level]->dof_indices
                = std::vector<types::global_dof_index> (next_free_dof,
                                                        numbers::invalid_dof_index);
              dof_handler.levels[level]->cell_dof_indices_cache
                = std::vector<types::global_dof_index> (cache_size,
                                                        numbers::invalid_dof_index);
            }

          // safety check: make sure that the number of DoFs we
          // allocated is actually correct (above we have also set the
          // dof_*_offsets field, so we couldn't use this simpler
          // algorithm)
#ifdef DEBUG
          for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
            {
              types::global_dof_index counter = 0;
              typename HpDoFHandler<dim,spacedim>::active_cell_iterator
              cell=dof_handler.begin_active(level),
              endc=dof_handler.end_active(level);
              for (; cell!=endc; ++cell)
                if (!cell->has_children()
                    &&
                    !cell->is_artificial())
                  counter += cell->get_fe().template n_dofs_per_object<dim>();

              Assert (dof_handler.levels[level]->dof_indices.size() == counter,
                      ExcInternalError());

              // also check that the number of unassigned slots in the dof_offsets
              // equals the number of cells on that level minus the number of
              // active, non-artificial cells (because these are exactly the cells
              // on which we do something)
              unsigned int n_active_non_artificial_cells = 0;
              for (cell=dof_handler.begin_active(level); cell!=endc; ++cell)
                if (!cell->has_children()
                    &&
                    !cell->is_artificial())
                  ++n_active_non_artificial_cells;

              Assert (static_cast<unsigned int>
                      (std::count (dof_handler.levels[level]->dof_offsets.begin(),
                                   dof_handler.levels[level]->dof_offsets.end(),
                                   (DoFLevel::offset_type)(-1)))
                      ==
                      dof_handler.tria->n_raw_cells(level) - n_active_non_artificial_cells,
                      ExcInternalError());
            }
#endif
        }




        /**
         * Do that part of reserving space that pertains to faces,
         * since this is the same in all space dimensions.
         */
        template <int dim, int spacedim>
        static
        void
        reserve_space_faces (DoFHandler<dim,spacedim> &dof_handler)
        {
          // make the code generic between lines and quads
          std::vector<unsigned int> &face_dof_offsets
            = (dim == 2
               ?
               reinterpret_cast<dealii::internal::hp::DoFIndicesOnFaces<2>&>(*dof_handler.faces).lines.dof_offsets
               :
               reinterpret_cast<dealii::internal::hp::DoFIndicesOnFaces<3>&>(*dof_handler.faces).quads.dof_offsets);

          std::vector<types::global_dof_index> &face_dof_indices
            = (dim == 2
               ?
               reinterpret_cast<dealii::internal::hp::DoFIndicesOnFaces<2>&>(*dof_handler.faces).lines.dofs
               :
               reinterpret_cast<dealii::internal::hp::DoFIndicesOnFaces<3>&>(*dof_handler.faces).quads.dofs);

          // FACE DOFS
          //
          // count face dofs, then allocate as much space
          // as we need and prime the linked list for faces (see the
          // description in hp::DoFLevel) with the indices we will
          // need. note that our task is more complicated than for the
          // cell case above since two adjacent cells may have different
          // active_fe_indices, in which case we need to allocate
          // *two* sets of face dofs for the same face
          //
          // the way we do things is that we loop over all active
          // cells (these are the ones that have DoFs only
          // anyway) and all their faces. We note in the
          // user flags whether we have previously visited a face and
          // if so skip it (consequently, we have to save and later
          // restore the face flags)
          {
            std::vector<bool> saved_face_user_flags;
            switch (dim)
              {
              case 2:
              {
                const_cast<dealii::Triangulation<dim,spacedim>&>(*dof_handler.tria)
                .save_user_flags_line (saved_face_user_flags);
                const_cast<dealii::Triangulation<dim,spacedim>&>(*dof_handler.tria)
                .clear_user_flags_line ();

                break;
              }

              case 3:
              {
                const_cast<dealii::Triangulation<dim,spacedim>&>(*dof_handler.tria)
                .save_user_flags_quad (saved_face_user_flags);
                const_cast<dealii::Triangulation<dim,spacedim>&>(*dof_handler.tria)
                .clear_user_flags_quad ();

                break;
              }

              default:
                Assert (false, ExcNotImplemented());
              }

            // an array to hold how many slots (see the hp::DoFLevel
            // class) we will have to store on each level
            unsigned int n_face_slots = 0;

            for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
                 cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
              if (! cell->is_artificial())
                for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
                  if (! cell->face(face)->user_flag_set())
                    {
                      // ok, face has not been visited. so we need to
                      // allocate space for it. let's see how much we
                      // need: we need one set if a) there is no
                      // neighbor behind this face, or b) the neighbor
                      // is either coarser or finer than we are, or c)
                      // the neighbor is artificial, or d) the neighbor
                      // is neither coarser nor finer, nor is artificial,
                      // and just so happens to have the same active_fe_index :
                      if (cell->at_boundary(face)
                          ||
                          cell->face(face)->has_children()
                          ||
                          cell->neighbor_is_coarser(face)
                          ||
                          (!cell->at_boundary(face)
                           &&
                           cell->neighbor(face)->is_artificial())
                          ||
                          (!cell->at_boundary(face)
                           &&
                           !cell->neighbor(face)->is_artificial()
                           &&
                           (cell->active_fe_index() == cell->neighbor(face)->active_fe_index())))
                        // ok, one set of dofs. that makes one active_fe_index, 1
                        // times dofs_per_face dofs, and one stop index
                        n_face_slots
                        += 1 + dof_handler.get_fe(cell->active_fe_index()).template n_dofs_per_object<dim-1>() + 1;

                      // otherwise we do indeed need two sets, i.e. two
                      // active_fe_indices, two sets of dofs, and one stop index:
                      else
                        n_face_slots
                        += (2
                            +
                            dof_handler.get_fe(cell->active_fe_index()).template n_dofs_per_object<dim-1>()
                            +
                            dof_handler.get_fe(cell->neighbor(face)->active_fe_index())
                            .template n_dofs_per_object<dim-1>()
                            +
                            1);

                      // mark this face as visited
                      cell->face(face)->set_user_flag ();
                    }

            // now that we know how many face dofs we will have to
            // have on each level, allocate the memory. note that we
            // allocate offsets for all faces, though only the active
            // ones will have a non-invalid value later on
            face_dof_offsets
              = std::vector<unsigned int> (dof_handler.tria->n_raw_faces(),
                                           (unsigned int)(-1));
            face_dof_indices
              = std::vector<types::global_dof_index> (n_face_slots,
                                                      numbers::invalid_dof_index);

            // with the memory now allocated, loop over the
            // dof_handler.cells again and prime the _offset values as
            // well as the fe_index fields
            switch (dim)
              {
              case 2:
              {
                const_cast<dealii::Triangulation<dim,spacedim>&>(*dof_handler.tria)
                .clear_user_flags_line ();

                break;
              }

              case 3:
              {
                const_cast<dealii::Triangulation<dim,spacedim>&>(*dof_handler.tria)
                .clear_user_flags_quad ();

                break;
              }

              default:
                Assert (false, ExcNotImplemented());
              }

            unsigned int next_free_face_slot = 0;

            for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
                 cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
              if (! cell->is_artificial())
                for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
                  if (! cell->face(face)->user_flag_set())
                    {
                      // same decision tree as before
                      if (cell->at_boundary(face)
                          ||
                          cell->face(face)->has_children()
                          ||
                          cell->neighbor_is_coarser(face)
                          ||
                          (!cell->at_boundary(face)
                           &&
                           cell->neighbor(face)->is_artificial())
                          ||
                          (!cell->at_boundary(face)
                           &&
                           !cell->neighbor(face)->is_artificial()
                           &&
                           (cell->active_fe_index() == cell->neighbor(face)->active_fe_index())))
                        {
                          face_dof_offsets[cell->face(face)->index()]
                            = next_free_face_slot;

                          // set first slot for this face to
                          // active_fe_index of this face
                          face_dof_indices[next_free_face_slot]
                            = cell->active_fe_index();

                          // the next dofs_per_face indices remain unset
                          // for the moment (i.e. at invalid_dof_index).
                          // following this comes the stop index, which
                          // also is invalid_dof_index and therefore
                          // does not have to be explicitly set

                          // finally, mark those slots as used
                          next_free_face_slot
                          += dof_handler.get_fe(cell->active_fe_index()).template n_dofs_per_object<dim-1>() + 2;
                        }
                      else
                        {
                          face_dof_offsets[cell->face(face)->index()]
                            = next_free_face_slot;

                          // set first slot for this face to
                          // active_fe_index of this face
                          face_dof_indices[next_free_face_slot]
                            = cell->active_fe_index();

                          // the next dofs_per_line/quad indices remain unset
                          // for the moment (i.e. at invalid_dof_index).
                          //
                          // then comes the fe_index for the neighboring
                          // cell:
                          face_dof_indices[next_free_face_slot
                                           +
                                           dof_handler.get_fe(cell->active_fe_index()).template n_dofs_per_object<dim-1>()
                                           +
                                           1]
                            = cell->neighbor(face)->active_fe_index();
                          // then again a set of dofs that we need not
                          // set right now
                          //
                          // following this comes the stop index, which
                          // also is invalid_dof_index and therefore
                          // does not have to be explicitly set

                          // finally, mark those slots as used
                          next_free_face_slot
                          += (dof_handler.get_fe(cell->active_fe_index()).template n_dofs_per_object<dim-1>()
                              +
                              dof_handler.get_fe(cell->neighbor(face)->active_fe_index())
                              .template n_dofs_per_object<dim-1>()
                              +
                              3);
                        }

                      // mark this face as visited
                      cell->face(face)->set_user_flag ();
                    }

            // we should have moved the cursor for each level to the
            // total number of dofs on that level. check that
            Assert (next_free_face_slot == n_face_slots,
                    ExcInternalError());

            // at the end, restore the user flags for the faces
            switch (dim)
              {
              case 2:
              {
                const_cast<dealii::Triangulation<dim,spacedim>&>(*dof_handler.tria)
                .load_user_flags_line (saved_face_user_flags);

                break;
              }

              case 3:
              {
                const_cast<dealii::Triangulation<dim,spacedim>&>(*dof_handler.tria)
                .load_user_flags_quad (saved_face_user_flags);

                break;
              }

              default:
                Assert (false, ExcNotImplemented());
              }

          }
        }


        /**
         * Reserve enough space in the <tt>levels[]</tt> objects to
         * store the numbers of the degrees of freedom needed for the
         * given element. The given element is that one which was
         * selected when calling @p distribute_dofs the last time.
         */
        template <int spacedim>
        static
        void
        reserve_space (DoFHandler<1,spacedim> &dof_handler)
        {
          const unsigned int dim = 1;

          typedef DoFHandler<dim,spacedim> BaseClass;

          Assert (dof_handler.finite_elements != nullptr,
                  typename BaseClass::ExcNoFESelected());
          Assert (dof_handler.finite_elements->size() > 0,
                  typename BaseClass::ExcNoFESelected());
          Assert (dof_handler.tria->n_levels() > 0,
                  ExcMessage("The current Triangulation must not be empty."));
          Assert (dof_handler.tria->n_levels() == dof_handler.levels.size (),
                  ExcInternalError ());

          reserve_space_release_space (dof_handler);

          reserve_space_cells (dof_handler);
          reserve_space_vertices (dof_handler);
        }



        template <int spacedim>
        static
        void
        reserve_space (DoFHandler<2,spacedim> &dof_handler)
        {
          const unsigned int dim = 2;

          typedef DoFHandler<dim,spacedim> BaseClass;

          Assert (dof_handler.finite_elements != nullptr,
                  typename BaseClass::ExcNoFESelected());
          Assert (dof_handler.finite_elements->size() > 0,
                  typename BaseClass::ExcNoFESelected());
          Assert (dof_handler.tria->n_levels() > 0,
                  ExcMessage("The current Triangulation must not be empty."));
          Assert (dof_handler.tria->n_levels() == dof_handler.levels.size (),
                  ExcInternalError ());

          reserve_space_release_space (dof_handler);

          // FIRST CELLS AND FACES
          reserve_space_cells (dof_handler);
          reserve_space_faces (dof_handler);

          // VERTEX DOFS
          reserve_space_vertices (dof_handler);
        }



        template <int spacedim>
        static
        void
        reserve_space (DoFHandler<3,spacedim> &dof_handler)
        {
          const unsigned int dim = 3;

          typedef DoFHandler<dim,spacedim> BaseClass;

          Assert (dof_handler.finite_elements != nullptr,
                  typename BaseClass::ExcNoFESelected());
          Assert (dof_handler.finite_elements->size() > 0,
                  typename BaseClass::ExcNoFESelected());
          Assert (dof_handler.tria->n_levels() > 0,
                  ExcMessage("The current Triangulation must not be empty."));
          Assert (dof_handler.tria->n_levels() == dof_handler.levels.size (),
                  ExcInternalError ());

          reserve_space_release_space (dof_handler);

          // FIRST CELLS AND FACES
          reserve_space_cells (dof_handler);
          reserve_space_faces (dof_handler);


          // LINE DOFS

          // the situation here is pretty much like with vertices:
          // there can be an arbitrary number of finite elements
          // associated with each line.
          //
          // the algorithm we use is somewhat similar to what we do in
          // reserve_space_vertices()
          if (true)
            {
              // what we do first is to set up an array in which we
              // record whether a line is associated with any of the
              // given fe's, by setting a bit. in a later step, we
              // then actually allocate memory for the required dofs
              std::vector<std::vector<bool> ><commit_msg>Remove a typedef thats unused when compiling in release mode.<commit_after>// ---------------------------------------------------------------------
//
// Copyright (C) 2003 - 2018 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.
//
// ---------------------------------------------------------------------

#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/dof_level.h>
#include <deal.II/hp/dof_faces.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_handler_policy.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_levels.h>
#include <deal.II/grid/tria.h>
#include <deal.II/fe/fe.h>
#include <deal.II/distributed/shared_tria.h>
#include <deal.II/distributed/tria.h>

#include <boost/serialization/array.hpp>

#include <set>
#include <algorithm>
#include <functional>

DEAL_II_NAMESPACE_OPEN

// The following is necessary for compilation under Visual Studio which is unable to correctly
// distinguish between dealii::DoFHandler and dealii::hp::DoFHandler.
// Plus it makes code in dof_handler.cc easier to read.
// Requires C++11 support which is in Visual Studio 2013 and newer.
#if defined(_MSC_VER) && (_MSC_VER >= 1800)
template <int dim, int spacedim> using HpDoFHandler = ::dealii::hp::DoFHandler<dim, spacedim>;
#else
// When using older Visual Studio or a different compiler just fall back.
#define HpDoFHandler DoFHandler
#endif

namespace parallel
{
  namespace distributed
  {
    template <int, int> class Triangulation;
  }
}




namespace internal
{
  namespace hp
  {
    namespace DoFHandler
    {
      // access class dealii::hp::DoFHandler instead of namespace
      // internal::hp::DoFHandler, etc
      using dealii::hp::DoFHandler;

      /**
       * A class with the same purpose as the similarly named class of the
       * Triangulation class. See there for more information.
       */
      struct Implementation
      {
        /**
         * Do that part of reserving space that pertains to releasing
         * the previously used memory.
         */
        template <int dim, int spacedim>
        static
        void
        reserve_space_release_space (DoFHandler<dim,spacedim> &dof_handler)
        {
          // Release all space except the active_fe_indices field
          // which we have to back up before
          {
            std::vector<std::vector<DoFLevel::active_fe_index_type> >
            active_fe_backup(dof_handler.levels.size ());
            for (unsigned int level = 0; level<dof_handler.levels.size (); ++level)
              active_fe_backup[level] = std::move(dof_handler.levels[level]->active_fe_indices);

            // delete all levels and set them up newly, since vectors
            // are troublesome if you want to change their size
            dof_handler.clear_space ();

            for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
              {
                dof_handler.levels.emplace_back (new internal::hp::DoFLevel);
                dof_handler.levels[level]->active_fe_indices = std::move(active_fe_backup[level]);
              }

            if (dim > 1)
              dof_handler.faces.reset (new internal::hp::DoFIndicesOnFaces<dim>);
          }
        }



        /**
         * Do that part of reserving space that pertains to vertices,
         * since this is the same in all space dimensions.
         */
        template <int dim, int spacedim>
        static
        void
        reserve_space_vertices (DoFHandler<dim,spacedim> &dof_handler)
        {
          // The final step in all of the reserve_space() functions is to set
          // up vertex dof information. since vertices are sequentially
          // numbered, what we do first is to set up an array in which
          // we record whether a vertex is associated with any of the
          // given fe's, by setting a bit. in a later step, we then
          // actually allocate memory for the required dofs
          //
          // in the following, we only need to consider vertices that are
          // adjacent to either a locally owned or a ghost cell; we never
          // store anything on vertices that are only surrounded by
          // artificial cells. so figure out that subset of vertices
          // first
          std::vector<bool> locally_used_vertices (dof_handler.tria->n_vertices(),
                                                   false);
          for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
               cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
            if (! cell->is_artificial())
              for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
                locally_used_vertices[cell->vertex_index(v)] = true;

          std::vector<std::vector<bool> >
          vertex_fe_association (dof_handler.finite_elements->size(),
                                 std::vector<bool> (dof_handler.tria->n_vertices(), false));

          for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
               cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
            if (! cell->is_artificial())
              for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
                vertex_fe_association[cell->active_fe_index()][cell->vertex_index(v)]
                  = true;

          // in debug mode, make sure that each vertex is associated
          // with at least one fe (note that except for unused
          // vertices, all vertices are actually active). this is of
          // course only true for vertices that are part of either
          // ghost or locally owned cells
#ifdef DEBUG
          for (unsigned int v=0; v<dof_handler.tria->n_vertices(); ++v)
            if (locally_used_vertices[v] == true)
              if (dof_handler.tria->vertex_used(v) == true)
                {
                  unsigned int fe=0;
                  for (; fe<dof_handler.finite_elements->size(); ++fe)
                    if (vertex_fe_association[fe][v] == true)
                      break;
                  Assert (fe != dof_handler.finite_elements->size(), ExcInternalError());
                }
#endif

          // next count how much memory we actually need. for each
          // vertex, we need one slot per fe to store the fe_index,
          // plus dofs_per_vertex for this fe. in addition, we need
          // one slot as the end marker for the fe_indices. at the
          // same time already fill the vertex_dof_offsets field
          dof_handler.vertex_dof_offsets.resize (dof_handler.tria->n_vertices(),
                                                 numbers::invalid_unsigned_int);

          unsigned int vertex_slots_needed = 0;
          for (unsigned int v=0; v<dof_handler.tria->n_vertices(); ++v)
            if (dof_handler.tria->vertex_used(v) == true)
              if (locally_used_vertices[v] == true)
                {
                  dof_handler.vertex_dof_offsets[v] = vertex_slots_needed;

                  for (unsigned int fe=0; fe<dof_handler.finite_elements->size(); ++fe)
                    if (vertex_fe_association[fe][v] == true)
                      vertex_slots_needed += dof_handler.get_fe(fe).dofs_per_vertex + 1;

                  // don't forget the end_marker:
                  ++vertex_slots_needed;
                }

          // now allocate the space we have determined we need, and
          // set up the linked lists for each of the vertices
          dof_handler.vertex_dofs.resize (vertex_slots_needed,
                                          numbers::invalid_dof_index);
          for (unsigned int v=0; v<dof_handler.tria->n_vertices(); ++v)
            if (dof_handler.tria->vertex_used(v) == true)
              if (locally_used_vertices[v] == true)
                {
                  unsigned int current_index = dof_handler.vertex_dof_offsets[v];
                  for (unsigned int fe=0; fe<dof_handler.finite_elements->size(); ++fe)
                    if (vertex_fe_association[fe][v] == true)
                      {
                        // if this vertex uses this fe, then set the
                        // fe_index and move the pointer ahead
                        dof_handler.vertex_dofs[current_index] = fe;
                        current_index += dof_handler.get_fe(fe).dofs_per_vertex + 1;
                      }
                  // finally place the end marker
                  dof_handler.vertex_dofs[current_index] = numbers::invalid_dof_index;
                }
        }


        /**
         * Do that part of reserving space that pertains to cells,
         * since this is the same in all space dimensions.
         */
        template <int dim, int spacedim>
        static
        void
        reserve_space_cells (DoFHandler<dim,spacedim> &dof_handler)
        {
          // count how much space we need on each level for the cell
          // dofs and set the dof_*_offsets data. initially set the
          // latter to an invalid index, and only later set it to
          // something reasonable for active dof_handler.cells
          //
          // note that for dof_handler.cells, the situation is simpler
          // than for other (lower dimensional) objects since exactly
          // one finite element is used for it
          for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
            {
              dof_handler.levels[level]->dof_offsets
                = std::vector<DoFLevel::offset_type> (dof_handler.tria->n_raw_cells(level),
                                                      (DoFLevel::offset_type)(-1));
              dof_handler.levels[level]->cell_cache_offsets
                = std::vector<DoFLevel::offset_type> (dof_handler.tria->n_raw_cells(level),
                                                      (DoFLevel::offset_type)(-1));

              types::global_dof_index next_free_dof = 0;
              types::global_dof_index cache_size = 0;
              typename HpDoFHandler<dim,spacedim>::active_cell_iterator
              cell=dof_handler.begin_active(level),
              endc=dof_handler.end_active(level);
              for (; cell!=endc; ++cell)
                if (!cell->has_children()
                    &&
                    !cell->is_artificial())
                  {
                    dof_handler.levels[level]->dof_offsets[cell->index()] = next_free_dof;
                    next_free_dof += cell->get_fe().template n_dofs_per_object<dim>();

                    dof_handler.levels[level]->cell_cache_offsets[cell->index()] = cache_size;
                    cache_size += cell->get_fe().dofs_per_cell;
                  }

              dof_handler.levels[level]->dof_indices
                = std::vector<types::global_dof_index> (next_free_dof,
                                                        numbers::invalid_dof_index);
              dof_handler.levels[level]->cell_dof_indices_cache
                = std::vector<types::global_dof_index> (cache_size,
                                                        numbers::invalid_dof_index);
            }

          // safety check: make sure that the number of DoFs we
          // allocated is actually correct (above we have also set the
          // dof_*_offsets field, so we couldn't use this simpler
          // algorithm)
#ifdef DEBUG
          for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
            {
              types::global_dof_index counter = 0;
              typename HpDoFHandler<dim,spacedim>::active_cell_iterator
              cell=dof_handler.begin_active(level),
              endc=dof_handler.end_active(level);
              for (; cell!=endc; ++cell)
                if (!cell->has_children()
                    &&
                    !cell->is_artificial())
                  counter += cell->get_fe().template n_dofs_per_object<dim>();

              Assert (dof_handler.levels[level]->dof_indices.size() == counter,
                      ExcInternalError());

              // also check that the number of unassigned slots in the dof_offsets
              // equals the number of cells on that level minus the number of
              // active, non-artificial cells (because these are exactly the cells
              // on which we do something)
              unsigned int n_active_non_artificial_cells = 0;
              for (cell=dof_handler.begin_active(level); cell!=endc; ++cell)
                if (!cell->has_children()
                    &&
                    !cell->is_artificial())
                  ++n_active_non_artificial_cells;

              Assert (static_cast<unsigned int>
                      (std::count (dof_handler.levels[level]->dof_offsets.begin(),
                                   dof_handler.levels[level]->dof_offsets.end(),
                                   (DoFLevel::offset_type)(-1)))
                      ==
                      dof_handler.tria->n_raw_cells(level) - n_active_non_artificial_cells,
                      ExcInternalError());
            }
#endif
        }




        /**
         * Do that part of reserving space that pertains to faces,
         * since this is the same in all space dimensions.
         */
        template <int dim, int spacedim>
        static
        void
        reserve_space_faces (DoFHandler<dim,spacedim> &dof_handler)
        {
          // make the code generic between lines and quads
          std::vector<unsigned int> &face_dof_offsets
            = (dim == 2
               ?
               reinterpret_cast<dealii::internal::hp::DoFIndicesOnFaces<2>&>(*dof_handler.faces).lines.dof_offsets
               :
               reinterpret_cast<dealii::internal::hp::DoFIndicesOnFaces<3>&>(*dof_handler.faces).quads.dof_offsets);

          std::vector<types::global_dof_index> &face_dof_indices
            = (dim == 2
               ?
               reinterpret_cast<dealii::internal::hp::DoFIndicesOnFaces<2>&>(*dof_handler.faces).lines.dofs
               :
               reinterpret_cast<dealii::internal::hp::DoFIndicesOnFaces<3>&>(*dof_handler.faces).quads.dofs);

          // FACE DOFS
          //
          // count face dofs, then allocate as much space
          // as we need and prime the linked list for faces (see the
          // description in hp::DoFLevel) with the indices we will
          // need. note that our task is more complicated than for the
          // cell case above since two adjacent cells may have different
          // active_fe_indices, in which case we need to allocate
          // *two* sets of face dofs for the same face
          //
          // the way we do things is that we loop over all active
          // cells (these are the ones that have DoFs only
          // anyway) and all their faces. We note in the
          // user flags whether we have previously visited a face and
          // if so skip it (consequently, we have to save and later
          // restore the face flags)
          {
            std::vector<bool> saved_face_user_flags;
            switch (dim)
              {
              case 2:
              {
                const_cast<dealii::Triangulation<dim,spacedim>&>(*dof_handler.tria)
                .save_user_flags_line (saved_face_user_flags);
                const_cast<dealii::Triangulation<dim,spacedim>&>(*dof_handler.tria)
                .clear_user_flags_line ();

                break;
              }

              case 3:
              {
                const_cast<dealii::Triangulation<dim,spacedim>&>(*dof_handler.tria)
                .save_user_flags_quad (saved_face_user_flags);
                const_cast<dealii::Triangulation<dim,spacedim>&>(*dof_handler.tria)
                .clear_user_flags_quad ();

                break;
              }

              default:
                Assert (false, ExcNotImplemented());
              }

            // an array to hold how many slots (see the hp::DoFLevel
            // class) we will have to store on each level
            unsigned int n_face_slots = 0;

            for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
                 cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
              if (! cell->is_artificial())
                for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
                  if (! cell->face(face)->user_flag_set())
                    {
                      // ok, face has not been visited. so we need to
                      // allocate space for it. let's see how much we
                      // need: we need one set if a) there is no
                      // neighbor behind this face, or b) the neighbor
                      // is either coarser or finer than we are, or c)
                      // the neighbor is artificial, or d) the neighbor
                      // is neither coarser nor finer, nor is artificial,
                      // and just so happens to have the same active_fe_index :
                      if (cell->at_boundary(face)
                          ||
                          cell->face(face)->has_children()
                          ||
                          cell->neighbor_is_coarser(face)
                          ||
                          (!cell->at_boundary(face)
                           &&
                           cell->neighbor(face)->is_artificial())
                          ||
                          (!cell->at_boundary(face)
                           &&
                           !cell->neighbor(face)->is_artificial()
                           &&
                           (cell->active_fe_index() == cell->neighbor(face)->active_fe_index())))
                        // ok, one set of dofs. that makes one active_fe_index, 1
                        // times dofs_per_face dofs, and one stop index
                        n_face_slots
                        += 1 + dof_handler.get_fe(cell->active_fe_index()).template n_dofs_per_object<dim-1>() + 1;

                      // otherwise we do indeed need two sets, i.e. two
                      // active_fe_indices, two sets of dofs, and one stop index:
                      else
                        n_face_slots
                        += (2
                            +
                            dof_handler.get_fe(cell->active_fe_index()).template n_dofs_per_object<dim-1>()
                            +
                            dof_handler.get_fe(cell->neighbor(face)->active_fe_index())
                            .template n_dofs_per_object<dim-1>()
                            +
                            1);

                      // mark this face as visited
                      cell->face(face)->set_user_flag ();
                    }

            // now that we know how many face dofs we will have to
            // have on each level, allocate the memory. note that we
            // allocate offsets for all faces, though only the active
            // ones will have a non-invalid value later on
            face_dof_offsets
              = std::vector<unsigned int> (dof_handler.tria->n_raw_faces(),
                                           (unsigned int)(-1));
            face_dof_indices
              = std::vector<types::global_dof_index> (n_face_slots,
                                                      numbers::invalid_dof_index);

            // with the memory now allocated, loop over the
            // dof_handler.cells again and prime the _offset values as
            // well as the fe_index fields
            switch (dim)
              {
              case 2:
              {
                const_cast<dealii::Triangulation<dim,spacedim>&>(*dof_handler.tria)
                .clear_user_flags_line ();

                break;
              }

              case 3:
              {
                const_cast<dealii::Triangulation<dim,spacedim>&>(*dof_handler.tria)
                .clear_user_flags_quad ();

                break;
              }

              default:
                Assert (false, ExcNotImplemented());
              }

            unsigned int next_free_face_slot = 0;

            for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
                 cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
              if (! cell->is_artificial())
                for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
                  if (! cell->face(face)->user_flag_set())
                    {
                      // same decision tree as before
                      if (cell->at_boundary(face)
                          ||
                          cell->face(face)->has_children()
                          ||
                          cell->neighbor_is_coarser(face)
                          ||
                          (!cell->at_boundary(face)
                           &&
                           cell->neighbor(face)->is_artificial())
                          ||
                          (!cell->at_boundary(face)
                           &&
                           !cell->neighbor(face)->is_artificial()
                           &&
                           (cell->active_fe_index() == cell->neighbor(face)->active_fe_index())))
                        {
                          face_dof_offsets[cell->face(face)->index()]
                            = next_free_face_slot;

                          // set first slot for this face to
                          // active_fe_index of this face
                          face_dof_indices[next_free_face_slot]
                            = cell->active_fe_index();

                          // the next dofs_per_face indices remain unset
                          // for the moment (i.e. at invalid_dof_index).
                          // following this comes the stop index, which
                          // also is invalid_dof_index and therefore
                          // does not have to be explicitly set

                          // finally, mark those slots as used
                          next_free_face_slot
                          += dof_handler.get_fe(cell->active_fe_index()).template n_dofs_per_object<dim-1>() + 2;
                        }
                      else
                        {
                          face_dof_offsets[cell->face(face)->index()]
                            = next_free_face_slot;

                          // set first slot for this face to
                          // active_fe_index of this face
                          face_dof_indices[next_free_face_slot]
                            = cell->active_fe_index();

                          // the next dofs_per_line/quad indices remain unset
                          // for the moment (i.e. at invalid_dof_index).
                          //
                          // then comes the fe_index for the neighboring
                          // cell:
                          face_dof_indices[next_free_face_slot
                                           +
                                           dof_handler.get_fe(cell->active_fe_index()).template n_dofs_per_object<dim-1>()
                                           +
                                           1]
                            = cell->neighbor(face)->active_fe_index();
                          // then again a set of dofs that we need not
                          // set right now
                          //
                          // following this comes the stop index, which
                          // also is invalid_dof_index and therefore
                          // does not have to be explicitly set

                          // finally, mark those slots as used
                          next_free_face_slot
                          += (dof_handler.get_fe(cell->active_fe_index()).template n_dofs_per_object<dim-1>()
                              +
                              dof_handler.get_fe(cell->neighbor(face)->active_fe_index())
                              .template n_dofs_per_object<dim-1>()
                              +
                              3);
                        }

                      // mark this face as visited
                      cell->face(face)->set_user_flag ();
                    }

            // we should have moved the cursor for each level to the
            // total number of dofs on that level. check that
            Assert (next_free_face_slot == n_face_slots,
                    ExcInternalError());

            // at the end, restore the user flags for the faces
            switch (dim)
              {
              case 2:
              {
                const_cast<dealii::Triangulation<dim,spacedim>&>(*dof_handler.tria)
                .load_user_flags_line (saved_face_user_flags);

                break;
              }

              case 3:
              {
                const_cast<dealii::Triangulation<dim,spacedim>&>(*dof_handler.tria)
                .load_user_flags_quad (saved_face_user_flags);

                break;
              }

              default:
                Assert (false, ExcNotImplemented());
              }

          }
        }


        /**
         * Reserve enough space in the <tt>levels[]</tt> objects to
         * store the numbers of the degrees of freedom needed for the
         * given element. The given element is that one which was
         * selected when calling @p distribute_dofs the last time.
         */
        template <int spacedim>
        static
        void
        reserve_space (DoFHandler<1,spacedim> &dof_handler)
        {
          const unsigned int dim = 1;

          Assert (dof_handler.finite_elements != nullptr,
                  (typename DoFHandler<dim,spacedim>::ExcNoFESelected()));
          Assert (dof_handler.finite_elements->size() > 0,
                  (typename DoFHandler<dim,spacedim>::ExcNoFESelected()));
          Assert (dof_handler.tria->n_levels() > 0,
                  ExcMessage("The current Triangulation must not be empty."));
          Assert (dof_handler.tria->n_levels() == dof_handler.levels.size (),
                  ExcInternalError ());

          reserve_space_release_space (dof_handler);

          reserve_space_cells (dof_handler);
          reserve_space_vertices (dof_handler);
        }



        template <int spacedim>
        static
        void
        reserve_space (DoFHandler<2,spacedim> &dof_handler)
        {
          const unsigned int dim = 2;

          Assert (dof_handler.finite_elements != nullptr,
                  (typename DoFHandler<dim,spacedim>::ExcNoFESelected()));
          Assert (dof_handler.finite_elements->size() > 0,
                  (typename DoFHandler<dim,spacedim>::ExcNoFESelected()));
          Assert (dof_handler.tria->n_levels() > 0,
                  ExcMessage("The current Triangulation must not be empty."));
          Assert (dof_handler.tria->n_levels() == dof_handler.levels.size (),
                  ExcInternalError ());

          reserve_space_release_space (dof_handler);

          // FIRST CELLS AND FACES
          reserve_space_cells (dof_handler);
          reserve_space_faces (dof_handler);

          // VERTEX DOFS
          reserve_space_vertices (dof_handler);
        }



        template <int spacedim>
        static
        void
        reserve_space (DoFHandler<3,spacedim> &dof_handler)
        {
          const unsigned int dim = 3;

          Assert (dof_handler.finite_elements != nullptr,
                  (typename DoFHandler<dim,spacedim>::ExcNoFESelected()));
          Assert (dof_handler.finite_elements->size() > 0,
                  (typename DoFHandler<dim,spacedim>::ExcNoFESelected()));
          Assert (dof_handler.tria->n_levels() > 0,
                  ExcMessage("The current Triangulation must not be empty."));
          Assert (dof_handler.tria->n_levels() == dof_handler.levels.size (),
                  ExcInternalError ());

          reserve_space_release_space (dof_handler);

          // FIRST CELLS AND FACES
          reserve_space_cells (dof_handler);
          reserve_space_faces (dof_handler);


          // LINE DOFS

          // the situation here is pretty much like with vertices:
          // there can be an arbitrary number of finite elements
          // associated with each line.
          //
          // the algorithm we use is somewhat similar to what we do in
          // reserve_space_vertices()
          if (true)
            {
              // what we do first is to set up an array in which we
              // record whether a line is associated with any of the
              // given fe's, by setting a bit. in a later step, we
              // then actually allocate memory for the required dofs
              std::vector<std::vector<bool> >
{'fraction_non_alphanumeric': 0.09406769626334219, 'fraction_numerical': 0.0029433157974863406, 'mean_word_length': 1.474902666722485, 'pattern_counts': {'":': 0, '<': 242, '<?xml version=': 0, '>': 450, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 70, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7427368', 'n_tokens_mistral': 15882, 'n_tokens_neox': 15133, 'n_words': 3904}