tria_objects.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_tria_objects_h
#define dealii_tria_objects_h

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

#include <deal.II/base/exceptions.h>
#include <deal.II/base/geometry_info.h>

#include <deal.II/grid/tria_object.h>

#include <vector>

DEAL_II_NAMESPACE_OPEN

// TODO: This should all be cleaned up. Currently, only a single
// function in the library makes use of the odd specializations, and
// this function is Triangulation::execute_refinement() in 3D. I
// assume, that the other refinement functions would profit from using
// next_free_single_object() and next_free_pair_object, but they seem
// to get around it.

// TODO: The TriaObjects class contains a std::vector<G>. This is only an
// efficient storage scheme if G is relatively well packed, i.e. it's not a
// bool and then an integer and then a double, etc. Verify that this is
// actually the case.

template <int dim, int spacedim>
class Triangulation;
template <class Accessor>
class TriaRawIterator;
template <int, int, int>
class TriaAccessor;

namespace internal
{
  namespace TriangulationImplementation
  {
    template <typename G>
    class TriaObjects
    {
    public:
      TriaObjects();

      std::vector<G> cells;

      std::vector<int> children;

      std::vector<RefinementCase<G::dimension>> refinement_cases;

      std::vector<bool> used;

      std::vector<bool> user_flags;


      struct BoundaryOrMaterialId
      {
        union
        {
          types::boundary_id boundary_id;
          types::material_id material_id;
        };


        BoundaryOrMaterialId();

        static std::size_t
        memory_consumption();

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

      std::vector<BoundaryOrMaterialId> boundary_or_material_id;

      std::vector<types::manifold_id> manifold_id;

      void
      reserve_space(const unsigned int new_objs_in_pairs,
                    const unsigned int new_objs_single = 0);

      template <int dim, int spacedim>
      ::TriaRawIterator<::TriaAccessor<G::dimension, dim, spacedim>>
      next_free_single_object(const ::Triangulation<dim, spacedim> &tria);

      template <int dim, int spacedim>
      ::TriaRawIterator<::TriaAccessor<G::dimension, dim, spacedim>>
      next_free_pair_object(const ::Triangulation<dim, spacedim> &tria);

      template <int dim, int spacedim>
      typename ::Triangulation<dim, spacedim>::raw_hex_iterator
      next_free_hex(const ::Triangulation<dim, spacedim> &tria,
                    const unsigned int                          level);

      void
      clear();

      bool
      face_orientation(const unsigned int cell, const unsigned int face) const;


      void *&
      user_pointer(const unsigned int i);

      const void *
      user_pointer(const unsigned int i) const;

      unsigned int &
      user_index(const unsigned int i);

      unsigned int
      user_index(const unsigned int i) const;

      void
      clear_user_data(const unsigned int i);

      void
      clear_user_data();

      void
      clear_user_flags();

      void
      monitor_memory(const unsigned int true_dimension) const;

      std::size_t
      memory_consumption() const;

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

      DeclException2(ExcMemoryInexact,
                     int,
                     int,
                     << "The containers have sizes " << arg1 << " and " << arg2
                     << ", which is not as expected.");

      DeclException0(ExcPointerIndexClash);

    protected:
      unsigned int next_free_single;

      unsigned int next_free_pair;

      bool reverse_order_next_free_single;

      struct UserData
      {
        union
        {
          void *p;
          unsigned int i;
        };

        UserData()
        {
          p = nullptr;
        }

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

      enum UserDataType
      {
        data_unknown,
        data_pointer,
        data_index
      };


      std::vector<UserData> user_data;

      mutable UserDataType user_data_type;
    };

    class TriaObjectsHex : public TriaObjects<TriaObject<3>>
    {
    public:
      bool
      face_orientation(const unsigned int cell, const unsigned int face) const;


      std::vector<bool> face_orientations;

      std::vector<bool> face_flips;

      std::vector<bool> face_rotations;

      void
      reserve_space(const unsigned int new_objs);

      void
      clear();

      void
      monitor_memory(const unsigned int true_dimension) const;

      std::size_t
      memory_consumption() const;

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


    class TriaObjectsQuad3D : public TriaObjects<TriaObject<2>>
    {
    public:
      bool
      face_orientation(const unsigned int cell, const unsigned int face) const;


      std::vector<bool> line_orientations;

      void
      reserve_space(const unsigned int new_quads_in_pairs,
                    const unsigned int new_quads_single = 0);

      void
      clear();

      void
      monitor_memory(const unsigned int true_dimension) const;

      std::size_t
      memory_consumption() const;

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

    //----------------------------------------------------------------------//


    template <typename G>
    inline TriaObjects<G>::BoundaryOrMaterialId::BoundaryOrMaterialId()
    {
      material_id = numbers::invalid_material_id;
    }



    template <typename G>
    std::size_t
    TriaObjects<G>::BoundaryOrMaterialId::memory_consumption()
    {
      return sizeof(BoundaryOrMaterialId);
    }



    template <typename G>
    template <class Archive>
    void
    TriaObjects<G>::BoundaryOrMaterialId::serialize(
      Archive &ar,
      const unsigned int /*version*/)
    {
      // serialize this
      // structure by
      // writing and
      // reading the larger
      // of the two values,
      // in order to make
      // sure we get all
      // bits
      if (sizeof(material_id) > sizeof(boundary_id))
        ar &material_id;
      else
        ar &boundary_id;
    }


    template <typename G>
    inline bool
    TriaObjects<G>::face_orientation(const unsigned int,
                                     const unsigned int) const
    {
      return true;
    }


    template <typename G>
    inline void *&
    TriaObjects<G>::user_pointer(const unsigned int i)
    {
      Assert(user_data_type == data_unknown || user_data_type == data_pointer,
             ExcPointerIndexClash());
      user_data_type = data_pointer;

      Assert(i < user_data.size(), ExcIndexRange(i, 0, user_data.size()));
      return user_data[i].p;
    }


    template <typename G>
    inline const void *
    TriaObjects<G>::user_pointer(const unsigned int i) const
    {
      Assert(user_data_type == data_unknown || user_data_type == data_pointer,
             ExcPointerIndexClash());
      user_data_type = data_pointer;

      Assert(i < user_data.size(), ExcIndexRange(i, 0, user_data.size()));
      return user_data[i].p;
    }


    template <typename G>
    inline unsigned int &
    TriaObjects<G>::user_index(const unsigned int i)
    {
      Assert(user_data_type == data_unknown || user_data_type == data_index,
             ExcPointerIndexClash());
      user_data_type = data_index;

      Assert(i < user_data.size(), ExcIndexRange(i, 0, user_data.size()));
      return user_data[i].i;
    }


    template <typename G>
    inline void
    TriaObjects<G>::clear_user_data(const unsigned int i)
    {
      Assert(i < user_data.size(), ExcIndexRange(i, 0, user_data.size()));
      user_data[i].i = 0;
    }


    template <typename G>
    inline TriaObjects<G>::TriaObjects()
      : next_free_single(numbers::invalid_unsigned_int)
      , next_free_pair(numbers::invalid_unsigned_int)
      , reverse_order_next_free_single(false)
      , user_data_type(data_unknown)
    {}


    template <typename G>
    inline unsigned int
    TriaObjects<G>::user_index(const unsigned int i) const
    {
      Assert(user_data_type == data_unknown || user_data_type == data_index,
             ExcPointerIndexClash());
      user_data_type = data_index;

      Assert(i < user_data.size(), ExcIndexRange(i, 0, user_data.size()));
      return user_data[i].i;
    }


    template <typename G>
    inline void
    TriaObjects<G>::clear_user_data()
    {
      user_data_type = data_unknown;
      for (unsigned int i = 0; i < user_data.size(); ++i)
        user_data[i].p = nullptr;
    }


    template <typename G>
    inline void
    TriaObjects<G>::clear_user_flags()
    {
      user_flags.assign(user_flags.size(), false);
    }


    template <typename G>
    template <class Archive>
    void
    TriaObjects<G>::UserData::serialize(Archive &ar, const unsigned int)
    {
      // serialize this as an integer
      ar &i;
    }



    template <typename G>
    template <class Archive>
    void
    TriaObjects<G>::serialize(Archive &ar, const unsigned int)
    {
      ar &cells &children;
      ar &       refinement_cases;
      ar &       used;
      ar &       user_flags;
      ar &       boundary_or_material_id;
      ar &       manifold_id;
      ar &next_free_single &next_free_pair &reverse_order_next_free_single;
      ar &user_data &user_data_type;
    }


    template <class Archive>
    void
    TriaObjectsHex::serialize(Archive &ar, const unsigned int version)
    {
      this->TriaObjects<TriaObject<3>>::serialize(ar, version);

      ar &face_orientations &face_flips &face_rotations;
    }


    template <class Archive>
    void
    TriaObjectsQuad3D::serialize(Archive &ar, const unsigned int version)
    {
      this->TriaObjects<TriaObject<2>>::serialize(ar, version);

      ar &line_orientations;
    }


    //----------------------------------------------------------------------//

    inline bool
    TriaObjectsHex::face_orientation(const unsigned int cell,
                                     const unsigned int face) const
    {
      Assert(cell < face_orientations.size() / GeometryInfo<3>::faces_per_cell,
             ExcIndexRange(0,
                           cell,
                           face_orientations.size() /
                             GeometryInfo<3>::faces_per_cell));
      Assert(face < GeometryInfo<3>::faces_per_cell,
             ExcIndexRange(0, face, GeometryInfo<3>::faces_per_cell));

      return face_orientations[cell * GeometryInfo<3>::faces_per_cell + face];
    }

    //----------------------------------------------------------------------//

    inline bool
    TriaObjectsQuad3D::face_orientation(const unsigned int cell,
                                        const unsigned int face) const
    {
      return line_orientations[cell * GeometryInfo<2>::faces_per_cell + face];
    }


    //----------------------------------------------------------------------//

    template <class G>
    template <int dim, int spacedim>
    ::TriaRawIterator<::TriaAccessor<G::dimension, dim, spacedim>>
    TriaObjects<G>::next_free_single_object(
      const ::Triangulation<dim, spacedim> &tria)
    {
      // TODO: Think of a way to ensure that we are using the correct
      // triangulation, i.e. the one containing *this.

      int pos = next_free_single, last = used.size() - 1;
      if (!reverse_order_next_free_single)
        {
          // first sweep forward, only use really single slots, do not use
          // pair slots
          for (; pos < last; ++pos)
            if (!used[pos])
              if (used[++pos])
                {
                  // this was a single slot
                  pos -= 1;
                  break;
                }
          if (pos >= last)
            {
              reverse_order_next_free_single = true;
              next_free_single               = used.size() - 1;
              pos                            = used.size() - 1;
            }
          else
            next_free_single = pos + 1;
        }

      if (reverse_order_next_free_single)
        {
          // second sweep, use all slots, even
          // in pairs
          for (; pos >= 0; --pos)
            if (!used[pos])
              break;
          if (pos > 0)
            next_free_single = pos - 1;
          else
            // no valid single object anymore
            return ::TriaRawIterator<
              ::TriaAccessor<G::dimension, dim, spacedim>>(&tria, -1, -1);
        }

      return ::TriaRawIterator<
        ::TriaAccessor<G::dimension, dim, spacedim>>(&tria, 0, pos);
    }



    template <class G>
    template <int dim, int spacedim>
    ::TriaRawIterator<::TriaAccessor<G::dimension, dim, spacedim>>
    TriaObjects<G>::next_free_pair_object(
      const ::Triangulation<dim, spacedim> &tria)
    {
      // TODO: Think of a way to ensure that we are using the correct
      // triangulation, i.e. the one containing *this.

      int pos = next_free_pair, last = used.size() - 1;
      for (; pos < last; ++pos)
        if (!used[pos])
          if (!used[++pos])
            {
              // this was a pair slot
              pos -= 1;
              break;
            }
      if (pos >= last)
        // no free slot
        return ::TriaRawIterator<
          ::TriaAccessor<G::dimension, dim, spacedim>>(&tria, -1, -1);
      else
        next_free_pair = pos + 2;

      return ::TriaRawIterator<
        ::TriaAccessor<G::dimension, dim, spacedim>>(&tria, 0, pos);
    }



    // declaration of explicit specializations

    template <>
    void
    TriaObjects<TriaObject<2>>::monitor_memory(const unsigned int) const;

  } // namespace TriangulationImplementation
} // namespace internal



DEAL_II_NAMESPACE_CLOSE

#endif