subdivision_levels.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/base_header.hpp>
#include <kernel/geometry/index_set.hpp>
#include <kernel/shape.hpp>
#include <kernel/util/assertion.hpp>
 
// includes, system
#include <functional>
 
namespace FEAT::Geometry
{
  template<int size_>
  struct SubdivisionLevelTuple
  {
    static_assert(size_ > 0, "invalid number of indices");
 
    static constexpr int size = size_;
 
    std::array<std::uint64_t, size_> levels;
 
    Index& operator[](Index idx)
    {
      ASSERT(idx < size);
      return levels[idx];
    }
 
    const Index& operator[](Index idx) const
    {
      ASSERT(idx < size);
      return levels[idx];
    }
 
    SubdivisionLevelTuple<size> decrement() const
    {
      SubdivisionLevelTuple<size> result;
      for(int i = 0; i < size; i++)
      {
        result[i] = levels[i] > 0 ? levels[i] - 1 : 0;
      }
      return result;
    }
  };
 
  class SubdivisionLevels
  {
    friend std::ostream& operator<<(std::ostream& /*unused*/, const SubdivisionLevels& /*unused*/);
 
  public:
    explicit SubdivisionLevels(Index num_vertices) : _levels(num_vertices)
    {
    }
 
    SubdivisionLevels(Index num_vertices, std::uint64_t default_level) : _levels(num_vertices, default_level)
    {
    }
 
    SubdivisionLevels(SubdivisionLevels&& other) noexcept :
      _levels(std::forward<std::vector<std::uint64_t>>(other._levels))
    {
    }
 
    SubdivisionLevels& operator=(SubdivisionLevels&& other) noexcept
    {
      // avoid self-move
      if(this == &other)
      {
        return *this;
      }
 
      _levels = std::forward<std::vector<std::uint64_t>>(other._levels);
 
      return *this;
    }
 
    SubdivisionLevels(const SubdivisionLevels& other) = delete;
    SubdivisionLevels& operator=(const SubdivisionLevels& other) = delete;
 
    ~SubdivisionLevels() = default;
 
    std::uint64_t& operator[](Index idx)
    {
      ASSERT(idx < _levels.size());
      return _levels[idx];
    }
 
    const std::uint64_t& operator[](Index idx) const
    {
      ASSERT(idx < _levels.size());
      return _levels[idx];
    }
 
    std::uint64_t maximum_level()
    {
      std::uint64_t result = 0;
      for(auto level : _levels)
      {
        result = std::max(result, level);
      }
      return result;
    }
 
    std::size_t size() const
    {
      return _levels.size();
    }
 
    template<typename MeshType>
    void set_per_element(const MeshType& mesh, std::function<std::uint64_t(Index)> f)
    {
      constexpr int cell_dim = MeshType::shape_dim;
 
      Index num_verts = mesh.get_num_entities(0);
      Index num_cells = mesh.get_num_entities(cell_dim);
 
      // Initialize levels to correct size
      _levels = std::vector<std::uint64_t>(num_verts, 0);
 
      for(Index cell(0); cell < num_cells; ++cell)
      {
        std::uint64_t cell_level = f(cell);
        set_level_for_element(mesh, cell, cell_level);
      }
    }
 
    template<typename MeshType>
    void set_level_for_element(const MeshType& mesh, Index element, std::uint64_t element_level)
    {
      constexpr int cell_dim = MeshType::shape_dim;
      constexpr int verts_per_cell = Shape::FaceTraits<typename MeshType::ShapeType, 0>::count;
 
      auto& v_at_c = mesh.template get_index_set<cell_dim, 0>();
      for(Index vert(0); vert < verts_per_cell; ++vert)
      {
        Index vert_idx = v_at_c(element, vert);
        _levels[vert_idx] = std::max(_levels[vert_idx], element_level);
      }
    }
 
    void merge(const SubdivisionLevels& other)
    {
      Index len = std::max(_levels.size(), other._levels.size());
      for(Index i(0); i < len; ++i)
      {
        _levels[i] += other._levels[i];
      }
    }
 
    template<int n>
    SubdivisionLevelTuple<n> get_levels(const IndexTuple<n>& indices) const
    {
      SubdivisionLevelTuple<n> result = {};
      for(int i = 0; i < indices.num_indices; i++)
      {
        result[i] = _levels[indices[i]];
      }
      return result;
    }
 
  private:
    std::vector<std::uint64_t> _levels;
  };
 
  inline std::ostream& operator<<(std::ostream& os, const SubdivisionLevels& sdls)
  {
    os << "SDLs {" << "\n";
    for(auto sdl : sdls._levels)
    {
      os << "\t" << sdl << ",\n";
    }
    os << "}\n";
    return os;
  }
} // namespace FEAT::Geometry