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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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