DOLFINx
DOLFINx 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/utils.h>
16 #include <dolfinx/graph/AdjacencyList.h>
17 #include <dolfinx/graph/partition.h>
18 #include <dolfinx/io/cells.h>
19 #include <map>
20 #include <memory>
21 #include <utility>
22 #include <vector>
23 #include <xtl/xspan.hpp>
24 
25 namespace dolfinx::mesh
26 {
27 
34 template <typename T>
35 class MeshTags
36 {
37 public:
45  template <typename U, typename V>
46  MeshTags(const std::shared_ptr<const Mesh>& mesh, int dim, U&& indices,
47  V&& values)
48  : _mesh(mesh), _dim(dim), _indices(std::forward<U>(indices)),
49  _values(std::forward<V>(values))
50  {
51  if (indices.size() != values.size())
52  {
53  throw std::runtime_error(
54  "Indices and values arrays must have same size.");
55  }
56 #ifdef DEBUG
57  if (!std::is_sorted(_indices.begin(), _indices.end()))
58  throw std::runtime_error("MeshTag data is not sorted");
59  if (std::adjacent_find(_indices.begin(), _indices.end()) != _indices.end())
60  throw std::runtime_error("MeshTag data has duplicates");
61 #endif
62  }
63 
65  MeshTags(const MeshTags& tags) = default;
66 
68  MeshTags(MeshTags&& tags) = default;
69 
71  ~MeshTags() = default;
72 
74  MeshTags& operator=(const MeshTags& tags) = default;
75 
77  MeshTags& operator=(MeshTags&& tags) = default;
78 
82  std::vector<std::int32_t> find(const T value) const
83  {
84  int n = std::count(_values.begin(), _values.end(), value);
85  std::vector<std::int32_t> indices(n);
86  int counter = 0;
87  for (std::int32_t i = 0; i < _values.size(); ++i)
88  {
89  if (_values[i] == value)
90  indices[counter++] = _indices[i];
91  }
92  return indices;
93  }
94 
97  const std::vector<std::int32_t>& indices() const { return _indices; }
98 
100  const std::vector<T>& values() const { return _values; }
101 
103  int dim() const { return _dim; }
104 
106  std::shared_ptr<const Mesh> mesh() const { return _mesh; }
107 
109  std::string name = "mesh_tags";
110 
112  std::size_t id() const { return _unique_id; }
113 
114 private:
115  // Unique identifier
116  std::size_t _unique_id = common::UniqueIdGenerator::id();
117 
118  // Associated mesh
119  std::shared_ptr<const Mesh> _mesh;
120 
121  // Topological dimension of tagged mesh entities
122  int _dim;
123 
124  // Local-to-process indices of tagged entities
125  std::vector<std::int32_t> _indices;
126 
127  // Values attached to entities
128  std::vector<T> _values;
129 };
130 
137 template <typename T>
139 create_meshtags(const std::shared_ptr<const mesh::Mesh>& mesh, const int dim,
140  const graph::AdjacencyList<std::int32_t>& entities,
141  const xtl::span<const T>& values)
142 {
143  assert(mesh);
144  if ((std::size_t)entities.num_nodes() != values.size())
145  throw std::runtime_error("Number of entities and values must match");
146 
147  // Tagged entity topological dimension
148  const auto map_e = mesh->topology().index_map(dim);
149  assert(map_e);
150 
151  auto e_to_v = mesh->topology().connectivity(dim, 0);
152  if (!e_to_v)
153  throw std::runtime_error("Missing entity-vertex connectivity.");
154 
155  const int num_vertices_per_entity = e_to_v->num_links(0);
156  const int num_entities_mesh = map_e->size_local() + map_e->num_ghosts();
157 
158  std::map<std::vector<std::int32_t>, std::int32_t> entity_key_to_index;
159  std::vector<std::int32_t> key(num_vertices_per_entity);
160  for (int e = 0; e < num_entities_mesh; ++e)
161  {
162  // Prepare a map from ordered local vertex indices to local entity
163  // number. This map is used to identify if received entity is owned
164  // or ghost
165  auto vertices = e_to_v->links(e);
166  std::copy(vertices.begin(), vertices.end(), key.begin());
167  std::sort(key.begin(), key.end());
168  entity_key_to_index.insert({key, e});
169  }
170 
171  // Iterate over all received entities. If entity is on this rank,
172  // store (local entity index, tag value)
173  std::vector<std::int32_t> indices_new;
174  std::vector<T> values_new;
175  std::vector<std::int32_t> entity(num_vertices_per_entity);
176  for (std::int32_t e = 0; e < entities.num_nodes(); ++e)
177  {
178  // This would fail for mixed cell type meshes
179  assert(num_vertices_per_entity == entities.num_links(e));
180  std::copy(entities.links(e).begin(), entities.links(e).end(),
181  entity.begin());
182  std::sort(entity.begin(), entity.end());
183 
184  if (const auto it = entity_key_to_index.find(entity);
185  it != entity_key_to_index.end())
186  {
187  indices_new.push_back(it->second);
188  values_new.push_back(values[e]);
189  }
190  }
191 
192  auto [indices_sorted, values_sorted]
193  = common::sort_unique(indices_new, values_new);
194  return mesh::MeshTags<T>(mesh, dim, std::move(indices_sorted),
195  std::move(values_sorted));
196 }
197 } // namespace dolfinx::mesh
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:47
xtl::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:130
std::int32_t num_nodes() const
Get the number of nodes.
Definition: AdjacencyList.h:115
int num_links(int node) const
Number of connections for given node.
Definition: AdjacencyList.h:120
A MeshTags are used to associate mesh entities with values. The entity index (local to process) ident...
Definition: MeshTags.h:36
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:103
std::shared_ptr< const Mesh > mesh() const
Return mesh.
Definition: MeshTags.h:106
std::size_t id() const
Unique ID of the object.
Definition: MeshTags.h:112
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:46
std::string name
Name.
Definition: MeshTags.h:109
~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:97
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:82
MeshTags & operator=(MeshTags &&tags)=default
Move assignment.
const std::vector< T > & values() const
Values attached to mesh entities.
Definition: MeshTags.h:100
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:29
Mesh data structures and algorithms on meshes.
Definition: DirichletBC.h:20
mesh::MeshTags< T > create_meshtags(const std::shared_ptr< const mesh::Mesh > &mesh, const int dim, const graph::AdjacencyList< std::int32_t > &entities, const xtl::span< const T > &values)
Create MeshTags from arrays.
Definition: MeshTags.h:139