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/utils.h>
16 #include <dolfinx/graph/AdjacencyList.h>
17 #include <dolfinx/graph/Partitioning.h>
18 #include <dolfinx/io/cells.h>
19 #include <map>
20 #include <memory>
21 #include <utility>
22 #include <vector>
23 
24 namespace dolfinx
25 {
26 namespace mesh
27 {
28 
35 template <typename T>
36 class MeshTags
37 {
38 public:
46  template <typename U, typename V>
47  MeshTags(const std::shared_ptr<const Mesh>& mesh, int dim, U&& indices,
48  V&& values)
49  : _mesh(mesh), _dim(dim), _indices(std::forward<U>(indices)),
50  _values(std::forward<V>(values))
51  {
52  if (indices.size() != values.size())
53  {
54  throw std::runtime_error(
55  "Indices and values arrays must have same size.");
56  }
57 #ifdef DEBUG
58  if (!std::is_sorted(_indices.begin(), _indices.end()))
59  throw std::runtime_error("MeshTag data is not sorted");
60  if (std::adjacent_find(_indices.begin(), _indices.end()) != _indices.end())
61  throw std::runtime_error("MeshTag data has duplicates");
62 #endif
63  }
64 
66  MeshTags(const MeshTags& tags) = default;
67 
69  MeshTags(MeshTags&& tags) = default;
70 
72  ~MeshTags() = default;
73 
75  MeshTags& operator=(const MeshTags& tags) = default;
76 
78  MeshTags& operator=(MeshTags&& tags) = default;
79 
83  const Eigen::Array<std::int32_t, Eigen::Dynamic, 1> find(const T value) const
84  {
85  int n = std::count(_values.begin(), _values.end(), value);
86  Eigen::Array<std::int32_t, Eigen::Dynamic, 1> indices(n);
87  int counter = 0;
88  for (std::int32_t i = 0; i < _values.size(); ++i)
89  {
90  if (_values[i] == value)
91  indices[counter++] = _indices[i];
92  }
93  return indices;
94  }
95 
98  const std::vector<std::int32_t>& indices() const { return _indices; }
99 
101  const std::vector<T>& values() const { return _values; }
102 
104  int dim() const { return _dim; }
105 
107  std::shared_ptr<const Mesh> mesh() const { return _mesh; }
108 
110  std::string name = "mesh_tags";
111 
113  std::size_t id() const { return _unique_id; }
114 
115 private:
116  // Unique identifier
117  std::size_t _unique_id = common::UniqueIdGenerator::id();
118 
119  // Associated mesh
120  std::shared_ptr<const Mesh> _mesh;
121 
122  // Topological dimension of tagged mesh entities
123  int _dim;
124 
125  // Local-to-process indices of tagged entities
126  std::vector<std::int32_t> _indices;
127 
128  // Values attached to entities
129  std::vector<T> _values;
130 };
131 
138 template <typename T>
140 create_meshtags(const std::shared_ptr<const mesh::Mesh>& mesh, const int dim,
141  const graph::AdjacencyList<std::int32_t>& entities,
142  const std::vector<T>& values)
143 {
144  assert(mesh);
145  if ((std::size_t)entities.num_nodes() != values.size())
146  throw std::runtime_error("Number of entities and values must match");
147 
148  // Tagged entity topological dimension
149  const auto map_e = mesh->topology().index_map(dim);
150  assert(map_e);
151 
152  auto e_to_v = mesh->topology().connectivity(dim, 0);
153  if (!e_to_v)
154  throw std::runtime_error("Missing entity-vertex connectivity.");
155 
156  const int num_vertices_per_entity = e_to_v->num_links(0);
157  const int num_entities_mesh = map_e->size_local() + map_e->num_ghosts();
158 
159  std::map<std::vector<std::int32_t>, std::int32_t> entity_key_to_index;
160  std::vector<std::int32_t> key(num_vertices_per_entity);
161  for (int e = 0; e < num_entities_mesh; ++e)
162  {
163  // Prepare a map from ordered local vertex indices to local entity number
164  // This map is used to identify if received entity is owned or ghost
165  auto vertices = e_to_v->links(e);
166  for (int v = 0; v < vertices.rows(); ++v)
167  key[v] = vertices[v];
168  std::sort(key.begin(), key.end());
169  entity_key_to_index.insert({key, e});
170  }
171 
172  // Iterate over all received entities. If entity is on this rank,
173  // store (local entity index, tag value)
174  std::vector<std::int32_t> indices_new;
175  std::vector<T> values_new;
176  std::vector<std::int32_t> entity(num_vertices_per_entity);
177 
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).data(),
183  entities.links(e).data() + num_vertices_per_entity,
184  entity.begin());
185  std::sort(entity.begin(), entity.end());
186 
187  if (const auto it = entity_key_to_index.find(entity);
188  it != entity_key_to_index.end())
189  {
190  indices_new.push_back(it->second);
191  values_new.push_back(values[e]);
192  }
193  }
194 
195  auto [indices_sorted, values_sorted]
196  = common::sort_unique(indices_new, values_new);
197  return mesh::MeshTags<T>(mesh, dim, std::move(indices_sorted),
198  std::move(values_sorted));
199 }
200 } // namespace mesh
201 } // namespace dolfinx
dolfinx::mesh::create_meshtags
mesh::MeshTags< T > create_meshtags(const std::shared_ptr< const mesh::Mesh > &mesh, const int dim, const graph::AdjacencyList< std::int32_t > &entities, const std::vector< T > &values)
Create MeshTags from arrays.
Definition: MeshTags.h:140
dolfinx::graph::AdjacencyList::num_links
int num_links(int node) const
Number of connections for given node.
Definition: AdjacencyList.h:143
dolfinx::graph::AdjacencyList::num_nodes
std::int32_t num_nodes() const
Number of nodes.
Definition: AdjacencyList.h:138
dolfinx::mesh::MeshTags::name
std::string name
Name.
Definition: MeshTags.h:110
dolfinx::graph::AdjacencyList::links
Eigen::Array< T, Eigen::Dynamic, 1 >::SegmentReturnType links(int node)
Links (edges) for given node.
Definition: AdjacencyList.h:153
dolfinx::mesh::MeshTags::dim
int dim() const
Return topological dimension of tagged entities.
Definition: MeshTags.h:104
dolfinx::mesh::MeshTags::mesh
std::shared_ptr< const Mesh > mesh() const
Return mesh.
Definition: MeshTags.h:107
dolfinx::common::sort_unique
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
dolfinx::graph::AdjacencyList
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition: AdjacencyList.h:28
dolfinx::mesh::MeshTags::MeshTags
MeshTags(MeshTags &&tags)=default
Move constructor.
dolfinx::mesh::MeshTags
A MeshTags are used to associate mesh entities with values. The entity index (local to process) ident...
Definition: MeshTags.h:37
dolfinx::mesh::MeshTags::operator=
MeshTags & operator=(const MeshTags &tags)=default
Move assignment.
dolfinx::mesh::MeshTags::id
std::size_t id() const
Unique ID of the object.
Definition: MeshTags.h:113
dolfinx::mesh::MeshTags::operator=
MeshTags & operator=(MeshTags &&tags)=default
Move assignment.
dolfinx::mesh::MeshTags::MeshTags
MeshTags(const MeshTags &tags)=default
Copy constructor.
dolfinx::mesh::MeshTags::indices
const std::vector< std::int32_t > & indices() const
Indices of tagged mesh entities (local-to-process). The indices are sorted.
Definition: MeshTags.h:98
dolfinx::mesh::MeshTags::values
const std::vector< T > & values() const
Values attached to mesh entities.
Definition: MeshTags.h:101
dolfinx::mesh::MeshTags::MeshTags
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:47
dolfinx::mesh::MeshTags::find
const Eigen::Array< std::int32_t, Eigen::Dynamic, 1 > find(const T value) const
Find all entities with a given tag value.
Definition: MeshTags.h:83
dolfinx::common::UniqueIdGenerator::id
static std::size_t id()
Generate a unique ID.
Definition: UniqueIdGenerator.cpp:22
dolfinx::mesh::MeshTags::~MeshTags
~MeshTags()=default
Destructor.