DOLFIN-X
DOLFIN-X C++ interface
MeshTags.h
1 // Copyright (C) 2020 Michal Habera
2 //
3 // This file is part of DOLFINX (https://www.fenicsproject.org)
4 //
5 // SPDX-License-Identifier: LGPL-3.0-or-later
6 
7 #pragma once
8 
9 #include "Geometry.h"
10 #include "Mesh.h"
11 #include "Topology.h"
12 #include <algorithm>
13 #include <dolfinx/common/IndexMap.h>
14 #include <dolfinx/common/UniqueIdGenerator.h>
15 #include <dolfinx/common/span.hpp>
16 #include <dolfinx/common/utils.h>
17 #include <dolfinx/graph/AdjacencyList.h>
18 #include <dolfinx/graph/partition.h>
19 #include <dolfinx/io/cells.h>
20 #include <map>
21 #include <memory>
22 #include <utility>
23 #include <vector>
24 
25 namespace dolfinx
26 {
27 namespace mesh
28 {
29 
36 template <typename T>
37 class MeshTags
38 {
39 public:
47  template <typename U, typename V>
48  MeshTags(const std::shared_ptr<const Mesh>& mesh, int dim, U&& indices,
49  V&& values)
50  : _mesh(mesh), _dim(dim), _indices(std::forward<U>(indices)),
51  _values(std::forward<V>(values))
52  {
53  if (indices.size() != values.size())
54  {
55  throw std::runtime_error(
56  "Indices and values arrays must have same size.");
57  }
58 #ifdef DEBUG
59  if (!std::is_sorted(_indices.begin(), _indices.end()))
60  throw std::runtime_error("MeshTag data is not sorted");
61  if (std::adjacent_find(_indices.begin(), _indices.end()) != _indices.end())
62  throw std::runtime_error("MeshTag data has duplicates");
63 #endif
64  }
65 
67  MeshTags(const MeshTags& tags) = default;
68 
70  MeshTags(MeshTags&& tags) = default;
71 
73  ~MeshTags() = default;
74 
76  MeshTags& operator=(const MeshTags& tags) = default;
77 
79  MeshTags& operator=(MeshTags&& tags) = default;
80 
84  std::vector<std::int32_t> find(const T value) const
85  {
86  int n = std::count(_values.begin(), _values.end(), value);
87  std::vector<std::int32_t> indices(n);
88  int counter = 0;
89  for (std::int32_t i = 0; i < _values.size(); ++i)
90  {
91  if (_values[i] == value)
92  indices[counter++] = _indices[i];
93  }
94  return indices;
95  }
96 
99  const std::vector<std::int32_t>& indices() const { return _indices; }
100 
102  const std::vector<T>& values() const { return _values; }
103 
105  int dim() const { return _dim; }
106 
108  std::shared_ptr<const Mesh> mesh() const { return _mesh; }
109 
111  std::string name = "mesh_tags";
112 
114  std::size_t id() const { return _unique_id; }
115 
116 private:
117  // Unique identifier
118  std::size_t _unique_id = common::UniqueIdGenerator::id();
119 
120  // Associated mesh
121  std::shared_ptr<const Mesh> _mesh;
122 
123  // Topological dimension of tagged mesh entities
124  int _dim;
125 
126  // Local-to-process indices of tagged entities
127  std::vector<std::int32_t> _indices;
128 
129  // Values attached to entities
130  std::vector<T> _values;
131 };
132 
139 template <typename T>
141 create_meshtags(const std::shared_ptr<const mesh::Mesh>& mesh, const int dim,
142  const graph::AdjacencyList<std::int32_t>& entities,
143  const tcb::span<const T>& values)
144 {
145  assert(mesh);
146  if ((std::size_t)entities.num_nodes() != values.size())
147  throw std::runtime_error("Number of entities and values must match");
148 
149  // Tagged entity topological dimension
150  const auto map_e = mesh->topology().index_map(dim);
151  assert(map_e);
152 
153  auto e_to_v = mesh->topology().connectivity(dim, 0);
154  if (!e_to_v)
155  throw std::runtime_error("Missing entity-vertex connectivity.");
156 
157  const int num_vertices_per_entity = e_to_v->num_links(0);
158  const int num_entities_mesh = map_e->size_local() + map_e->num_ghosts();
159 
160  std::map<std::vector<std::int32_t>, std::int32_t> entity_key_to_index;
161  std::vector<std::int32_t> key(num_vertices_per_entity);
162  for (int e = 0; e < num_entities_mesh; ++e)
163  {
164  // Prepare a map from ordered local vertex indices to local entity
165  // number. This map is used to identify if received entity is owned
166  // or ghost
167  auto vertices = e_to_v->links(e);
168  std::copy(vertices.begin(), vertices.end(), key.begin());
169  std::sort(key.begin(), key.end());
170  entity_key_to_index.insert({key, e});
171  }
172 
173  // Iterate over all received entities. If entity is on this rank,
174  // store (local entity index, tag value)
175  std::vector<std::int32_t> indices_new;
176  std::vector<T> values_new;
177  std::vector<std::int32_t> entity(num_vertices_per_entity);
178  for (Eigen::Index e = 0; e < entities.num_nodes(); ++e)
179  {
180  // This would fail for mixed cell type meshes
181  assert(num_vertices_per_entity == entities.num_links(e));
182  std::copy(entities.links(e).begin(), entities.links(e).end(),
183  entity.begin());
184  std::sort(entity.begin(), entity.end());
185 
186  if (const auto it = entity_key_to_index.find(entity);
187  it != entity_key_to_index.end())
188  {
189  indices_new.push_back(it->second);
190  values_new.push_back(values[e]);
191  }
192  }
193 
194  auto [indices_sorted, values_sorted]
195  = common::sort_unique(indices_new, values_new);
196  return mesh::MeshTags<T>(mesh, dim, std::move(indices_sorted),
197  std::move(values_sorted));
198 }
199 } // namespace mesh
200 } // namespace dolfinx
static std::size_t id()
Generate a unique ID.
Definition: UniqueIdGenerator.cpp:22
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition: AdjacencyList.h:46
tcb::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:128
std::int32_t num_nodes() const
Get the number of nodes.
Definition: AdjacencyList.h:113
int num_links(int node) const
Number of connections for given node.
Definition: AdjacencyList.h:118
A MeshTags are used to associate mesh entities with values. The entity index (local to process) ident...
Definition: MeshTags.h:38
MeshTags(const MeshTags &tags)=default
Copy constructor.
MeshTags & operator=(const MeshTags &tags)=default
Move assignment.
int dim() const
Return topological dimension of tagged entities.
Definition: MeshTags.h:105
std::shared_ptr< const Mesh > mesh() const
Return mesh.
Definition: MeshTags.h:108
std::size_t id() const
Unique ID of the object.
Definition: MeshTags.h:114
MeshTags(const std::shared_ptr< const Mesh > &mesh, int dim, U &&indices, V &&values)
Create from entities of given dimension on a mesh.
Definition: MeshTags.h:48
std::string name
Name.
Definition: MeshTags.h:111
~MeshTags()=default
Destructor.
const std::vector< std::int32_t > & indices() const
Indices of tagged mesh entities (local-to-process). The indices are sorted.
Definition: MeshTags.h:99
MeshTags(MeshTags &&tags)=default
Move constructor.
std::vector< std::int32_t > find(const T value) const
Find all entities with a given tag value.
Definition: MeshTags.h:84
MeshTags & operator=(MeshTags &&tags)=default
Move assignment.
const std::vector< T > & values() const
Values attached to mesh entities.
Definition: MeshTags.h:102
std::pair< U, V > sort_unique(const U &indices, const V &values)
Sort two arrays based on the values in array indices. Any duplicate indices and the corresponding val...
Definition: utils.h:30
mesh::MeshTags< T > create_meshtags(const std::shared_ptr< const mesh::Mesh > &mesh, const int dim, const graph::AdjacencyList< std::int32_t > &entities, const tcb::span< const T > &values)
Create MeshTags from arrays.
Definition: MeshTags.h:141