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

#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;

  }
}



DEAL_II_NAMESPACE_CLOSE

#endif