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