MeshTags.h

// Copyright (C) 2020-2022 Michal Habera and Garth N. Wells
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier:    LGPL-3.0-or-later
 
#pragma once
 
#include "Topology.h"
#include <algorithm>
#include <dolfinx/common/IndexMap.h>
#include <dolfinx/common/log.h>
#include <dolfinx/common/utils.h>
#include <dolfinx/graph/AdjacencyList.h>
#include <dolfinx/io/cells.h>
#include <memory>
#include <span>
#include <utility>
#include <vector>
 
namespace dolfinx::mesh
{
 
template <typename T>
class MeshTags
{
public:
  template <typename U, typename V>
    requires std::is_convertible_v<std::remove_cvref_t<U>,
                                   std::vector<std::int32_t>>
                 and std::is_convertible_v<std::remove_cvref_t<V>,
                                           std::vector<T>>
  MeshTags(std::shared_ptr<const Topology> topology, int dim, U&& indices,
           V&& values)
      : _topology(topology), _dim(dim), _indices(std::forward<U>(indices)),
        _values(std::forward<V>(values))
  {
    if (_indices.size() != _values.size())
    {
      throw std::runtime_error(
          "Indices and values arrays must have same size.");
    }
#ifndef NDEBUG
    if (!std::is_sorted(_indices.begin(), _indices.end()))
      throw std::runtime_error("MeshTag data is not sorted");
    if (std::adjacent_find(_indices.begin(), _indices.end()) != _indices.end())
      throw std::runtime_error("MeshTag data has duplicates");
#endif
  }
 
  MeshTags(const MeshTags& tags) = default;
 
  MeshTags(MeshTags&& tags) = default;
 
  ~MeshTags() = default;
 
  MeshTags& operator=(const MeshTags& tags) = default;
 
  MeshTags& operator=(MeshTags&& tags) = default;
 
  std::vector<std::int32_t> find(const T value) const
  {
    std::size_t n = std::count(_values.begin(), _values.end(), value);
    std::vector<std::int32_t> indices;
    indices.reserve(n);
    for (std::int32_t i = 0; i < _values.size(); ++i)
    {
      if (_values[i] == value)
        indices.push_back(_indices[i]);
    }
    return indices;
  }
 
  std::span<const std::int32_t> indices() const { return _indices; }
 
  std::span<const T> values() const { return _values; }
 
  int dim() const { return _dim; }
 
  std::shared_ptr<const Topology> topology() const { return _topology; }
 
  std::string name = "mesh_tags";
 
private:
  // Associated topology
  std::shared_ptr<const Topology> _topology;
 
  // Topological dimension of tagged topology entities
  int _dim;
 
  // Local-to-process indices of tagged entities
  std::vector<std::int32_t> _indices;
 
  // Values attached to entities
  std::vector<T> _values;
};
 
template <typename T>
MeshTags<T> create_meshtags(std::shared_ptr<const Topology> topology, int dim,
                            const graph::AdjacencyList<std::int32_t>& entities,
                            std::span<const T> values)
{
  LOG(INFO)
      << "Building MeshTags object from tagged entities (defined by vertices).";
 
  // Compute the indices of the topology entities (index is set to -1 if
  // it can't be found)
  assert(topology);
  const std::vector<std::int32_t> indices
      = entities_to_index(*topology, dim, entities.array());
  if (indices.size() != values.size())
  {
    throw std::runtime_error(
        "Duplicate mesh entities when building MeshTags object.");
  }
 
  // Sort the indices and values by indices
  auto [indices_sorted, values_sorted] = common::sort_unique(indices, values);
 
  // Remove any entities that were not found (these have an index of -1)
  auto it0 = std::lower_bound(indices_sorted.begin(), indices_sorted.end(), 0);
  std::size_t pos0 = std::distance(indices_sorted.begin(), it0);
  indices_sorted.erase(indices_sorted.begin(), it0);
  values_sorted.erase(values_sorted.begin(),
                      std::next(values_sorted.begin(), pos0));
 
  return MeshTags<T>(topology, dim, std::move(indices_sorted),
                     std::move(values_sorted));
}
} // namespace dolfinx::mesh