dof_faces.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2006 - 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_hp_dof_faces_h
#define dealii_hp_dof_faces_h


#include <deal.II/base/config.h>

#include <deal.II/hp/fe_collection.h>

#include <vector>

DEAL_II_NAMESPACE_OPEN

namespace internal
{
  namespace hp
  {
    template <int structdim>
    class DoFIndicesOnFacesOrEdges
    {
    public:
      std::vector<unsigned int> dof_offsets;

      std::vector<types::global_dof_index> dofs;

      template <int dim, int spacedim>
      void
      set_dof_index(const ::hp::DoFHandler<dim, spacedim> &dof_handler,
                    const unsigned int                           obj_index,
                    const unsigned int                           fe_index,
                    const unsigned int                           local_index,
                    const types::global_dof_index                global_index,
                    const unsigned int                           obj_level);

      template <int dim, int spacedim>
      types::global_dof_index
      get_dof_index(const ::hp::DoFHandler<dim, spacedim> &dof_handler,
                    const unsigned int                           obj_index,
                    const unsigned int                           fe_index,
                    const unsigned int                           local_index,
                    const unsigned int obj_level) const;

      template <int dim, int spacedim>
      unsigned int
      n_active_fe_indices(
        const ::hp::DoFHandler<dim, spacedim> &dof_handler,
        const unsigned int                           obj_index) const;

      template <int dim, int spacedim>
      types::global_dof_index
      nth_active_fe_index(
        const ::hp::DoFHandler<dim, spacedim> &dof_handler,
        const unsigned int                           obj_level,
        const unsigned int                           obj_index,
        const unsigned int                           n) const;

      template <int dim, int spacedim>
      bool
      fe_index_is_active(
        const ::hp::DoFHandler<dim, spacedim> &dof_handler,
        const unsigned int                           obj_index,
        const unsigned int                           fe_index,
        const unsigned int                           obj_level) const;

      std::size_t
      memory_consumption() const;

      template <class Archive>
      void
      serialize(Archive &ar, const unsigned int version);
    };



    template <int dim>
    class DoFIndicesOnFaces;


    template <>
    class DoFIndicesOnFaces<1>
    {
    public:
      std::size_t
      memory_consumption() const;

      template <class Archive>
      void
      serialize(Archive &ar, const unsigned int version);
    };

    template <>
    class DoFIndicesOnFaces<2>
    {
    public:
      internal::hp::DoFIndicesOnFacesOrEdges<1> lines;

      std::size_t
      memory_consumption() const;

      template <class Archive>
      void
      serialize(Archive &ar, const unsigned int version);
    };

    template <>
    class DoFIndicesOnFaces<3>
    {
    public:
      internal::hp::DoFIndicesOnFacesOrEdges<1> lines;

      internal::hp::DoFIndicesOnFacesOrEdges<2> quads;

      std::size_t
      memory_consumption() const;

      template <class Archive>
      void
      serialize(Archive &ar, const unsigned int version);
    };


    // --------------------- inline and template functions ------------------
    template <class Archive>
    void
    DoFIndicesOnFaces<1>::serialize(Archive &, const unsigned int)
    {}


    template <class Archive>
    void
    DoFIndicesOnFaces<2>::serialize(Archive &ar, const unsigned int)
    {
      ar &lines;
    }


    template <class Archive>
    void
    DoFIndicesOnFaces<3>::serialize(Archive &ar, const unsigned int)
    {
      ar &lines &quads;
    }

    template <int structdim>
    template <int dim, int spacedim>
    inline types::global_dof_index
    DoFIndicesOnFacesOrEdges<structdim>::get_dof_index(
      const ::hp::DoFHandler<dim, spacedim> &dof_handler,
      const unsigned int                           obj_index,
      const unsigned int                           fe_index,
      const unsigned int                           local_index,
      const unsigned int /*obj_level*/) const
    {
      Assert((fe_index !=
              ::hp::DoFHandler<dim, spacedim>::default_fe_index),
             ExcMessage("You need to specify a FE index when working "
                        "with hp DoFHandlers"));
      Assert(fe_index < dof_handler.get_fe_collection().size(),
             ExcIndexRange(fe_index,
                           0,
                           dof_handler.get_fe_collection().size()));
      Assert(local_index < dof_handler.get_fe(fe_index)
                             .template n_dofs_per_object<structdim>(),
             ExcIndexRange(local_index,
                           0,
                           dof_handler.get_fe(fe_index)
                             .template n_dofs_per_object<structdim>()));
      Assert(obj_index < dof_offsets.size(),
             ExcIndexRange(obj_index, 0, dof_offsets.size()));

      // make sure we are on an
      // object for which DoFs have
      // been allocated at all
      Assert(dof_offsets[obj_index] != numbers::invalid_unsigned_int,
             ExcMessage("You are trying to access degree of freedom "
                        "information for an object on which no such "
                        "information is available"));

      Assert(structdim < dim,
             ExcMessage("This object can not be used for cells."));

      // there may be multiple finite elements associated with
      // this object. hop along the list of index sets until we
      // find the one with the correct fe_index, and then poke
      // into that part. trigger an exception if we can't find a
      // set for this particular fe_index
      const types::global_dof_index  starting_offset = dof_offsets[obj_index];
      const types::global_dof_index *pointer         = &dofs[starting_offset];
      while (true)
        {
          Assert(*pointer != numbers::invalid_dof_index, ExcInternalError());
          if (*pointer == fe_index)
            return *(pointer + 1 + local_index);
          else
            pointer += static_cast<types::global_dof_index>(
              dof_handler.get_fe(*pointer)
                .template n_dofs_per_object<structdim>() +
              1);
        }
    }



    template <int structdim>
    template <int dim, int spacedim>
    inline void
    DoFIndicesOnFacesOrEdges<structdim>::set_dof_index(
      const ::hp::DoFHandler<dim, spacedim> &dof_handler,
      const unsigned int                           obj_index,
      const unsigned int                           fe_index,
      const unsigned int                           local_index,
      const types::global_dof_index                global_index,
      const unsigned int /*obj_level*/)
    {
      Assert((fe_index !=
              ::hp::DoFHandler<dim, spacedim>::default_fe_index),
             ExcMessage("You need to specify a FE index when working "
                        "with hp DoFHandlers"));
      Assert(fe_index < dof_handler.get_fe_collection().size(),
             ExcIndexRange(fe_index,
                           0,
                           dof_handler.get_fe_collection().size()));
      Assert(local_index < dof_handler.get_fe(fe_index)
                             .template n_dofs_per_object<structdim>(),
             ExcIndexRange(local_index,
                           0,
                           dof_handler.get_fe(fe_index)
                             .template n_dofs_per_object<structdim>()));
      Assert(obj_index < dof_offsets.size(),
             ExcIndexRange(obj_index, 0, dof_offsets.size()));

      // make sure we are on an
      // object for which DoFs have
      // been allocated at all
      Assert(dof_offsets[obj_index] != numbers::invalid_unsigned_int,
             ExcMessage("You are trying to access degree of freedom "
                        "information for an object on which no such "
                        "information is available"));

      Assert(structdim < dim,
             ExcMessage("This object can not be used for cells."));

      // there may be multiple finite elements associated with
      // this object.  hop along the list of index sets until we
      // find the one with the correct fe_index, and then poke
      // into that part. trigger an exception if we can't find a
      // set for this particular fe_index
      const types::global_dof_index starting_offset = dof_offsets[obj_index];
      types::global_dof_index *     pointer         = &dofs[starting_offset];
      while (true)
        {
          Assert(*pointer != numbers::invalid_dof_index, ExcInternalError());
          if (*pointer == fe_index)
            {
              *(pointer + 1 + local_index) = global_index;
              return;
            }
          else
            pointer += dof_handler.get_fe(*pointer)
                         .template n_dofs_per_object<structdim>() +
                       1;
        }
    }



    template <int structdim>
    template <int dim, int spacedim>
    inline unsigned int
    DoFIndicesOnFacesOrEdges<structdim>::n_active_fe_indices(
      const ::hp::DoFHandler<dim, spacedim> &dof_handler,
      const unsigned int                           obj_index) const
    {
      Assert(obj_index < dof_offsets.size(),
             ExcIndexRange(obj_index, 0, dof_offsets.size()));

      // make sure we are on an
      // object for which DoFs have
      // been allocated at all
      if (dof_offsets[obj_index] == numbers::invalid_unsigned_int)
        return 0;

      Assert(structdim < dim,
             ExcMessage("This object can not be used for cells."));

      // there may be multiple finite elements associated with this
      // object. hop along the list of index sets until we find the
      // one with the correct fe_index, and then poke into that
      // part. trigger an exception if we can't find a set for this
      // particular fe_index
      const unsigned int             starting_offset = dof_offsets[obj_index];
      const types::global_dof_index *pointer         = &dofs[starting_offset];
      unsigned int                   counter         = 0;
      while (true)
        {
          if (*pointer == numbers::invalid_dof_index)
            // end of list reached
            return counter;
          else
            {
              ++counter;
              pointer += dof_handler.get_fe(*pointer)
                           .template n_dofs_per_object<structdim>() +
                         1;
            }
        }
    }



    template <int structdim>
    template <int dim, int spacedim>
    inline types::global_dof_index
    DoFIndicesOnFacesOrEdges<structdim>::nth_active_fe_index(
      const ::hp::DoFHandler<dim, spacedim> &dof_handler,
      const unsigned int /*obj_level*/,
      const unsigned int obj_index,
      const unsigned int n) const
    {
      Assert(obj_index < dof_offsets.size(),
             ExcIndexRange(obj_index, 0, dof_offsets.size()));

      // make sure we are on an
      // object for which DoFs have
      // been allocated at all
      Assert(dof_offsets[obj_index] != numbers::invalid_unsigned_int,
             ExcMessage("You are trying to access degree of freedom "
                        "information for an object on which no such "
                        "information is available"));

      Assert(structdim < dim,
             ExcMessage("This object can not be used for cells."));

      Assert(n < n_active_fe_indices(dof_handler, obj_index),
             ExcIndexRange(n, 0, n_active_fe_indices(dof_handler, obj_index)));

      // there may be multiple finite elements associated with
      // this object. hop along the list of index sets until we
      // find the one with the correct fe_index, and then poke
      // into that part. trigger an exception if we can't find a
      // set for this particular fe_index
      const unsigned int             starting_offset = dof_offsets[obj_index];
      const types::global_dof_index *pointer         = &dofs[starting_offset];
      unsigned int                   counter         = 0;
      while (true)
        {
          Assert(*pointer != numbers::invalid_dof_index, ExcInternalError());

          const unsigned int fe_index = *pointer;

          Assert(fe_index < dof_handler.get_fe_collection().size(),
                 ExcInternalError());

          if (counter == n)
            return fe_index;

          ++counter;
          pointer += dof_handler.get_fe(fe_index)
                       .template n_dofs_per_object<structdim>() +
                     1;
        }
    }



    template <int structdim>
    template <int dim, int spacedim>
    inline bool
    DoFIndicesOnFacesOrEdges<structdim>::fe_index_is_active(
      const ::hp::DoFHandler<dim, spacedim> &dof_handler,
      const unsigned int                           obj_index,
      const unsigned int                           fe_index,
      const unsigned int /*obj_level*/) const
    {
      Assert(obj_index < dof_offsets.size(),
             ExcIndexRange(obj_index,
                           0,
                           static_cast<unsigned int>(dof_offsets.size())));
      Assert((fe_index !=
              ::hp::DoFHandler<dim, spacedim>::default_fe_index),
             ExcMessage("You need to specify a FE index when working "
                        "with hp DoFHandlers"));
      Assert(fe_index < dof_handler.get_fe_collection().size(),
             ExcIndexRange(fe_index,
                           0,
                           dof_handler.get_fe_collection().size()));

      // make sure we are on an
      // object for which DoFs have
      // been allocated at all
      Assert(dof_offsets[obj_index] != numbers::invalid_unsigned_int,
             ExcMessage("You are trying to access degree of freedom "
                        "information for an object on which no such "
                        "information is available"));

      Assert(structdim < dim,
             ExcMessage("This object can not be used for cells."));

      // there may be multiple finite elements associated with
      // this object. hop along the list of index sets until we
      // find the one with the correct fe_index, and then poke
      // into that part. trigger an exception if we can't find a
      // set for this particular fe_index
      const types::global_dof_index  starting_offset = dof_offsets[obj_index];
      const types::global_dof_index *pointer         = &dofs[starting_offset];
      while (true)
        {
          if (*pointer == numbers::invalid_dof_index)
            // end of list reached
            return false;
          else if (*pointer == fe_index)
            return true;
          else
            pointer += static_cast<types::global_dof_index>(
              dof_handler.get_fe(*pointer)
                .template n_dofs_per_object<structdim>() +
              1);
        }
    }

    template <int structdim>
    template <class Archive>
    void
    DoFIndicesOnFacesOrEdges<structdim>::serialize(Archive &ar,
                                                   const unsigned int)
    {
      ar &dofs;
      ar &dof_offsets;
    }


  } // namespace hp

} // namespace internal

DEAL_II_NAMESPACE_CLOSE

#endif