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 
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 
145 template <typename T>
147 create_meshtags(MPI_Comm comm, const std::shared_ptr<const mesh::Mesh>& mesh,
148  const mesh::CellType& tag_cell_type,
149  const Eigen::Array<std::int64_t, Eigen::Dynamic, Eigen::Dynamic,
150  Eigen::RowMajor>& entities,
151  const std::vector<T>& values)
152 {
153  // NOTE: Not yet working for higher-order geometries
154  //
155  // TODO: Avoid expensive-to-create std::vector<std::vector>>. Build
156  // AdjacencyList instead.
157 
158  assert(mesh);
159  if ((std::size_t)entities.rows() != values.size())
160  throw std::runtime_error("Number of entities and values must match");
161 
162  // Tagged entity topological dimension
163  const int e_dim = mesh::cell_dim(tag_cell_type);
164 
165  // -------------------
166  // 1. Send this rank's global "input" nodes indices to the
167  // 'postmaster' rank, and receive global "input" nodes for which
168  // this rank is the postmaster
169 
170  // Get "input" global node indices (as in the input file before any
171  // internal re-ordering)
172  const std::vector<std::int64_t>& nodes_g
173  = mesh->geometry().input_global_indices();
174 
175  // Send input global indices to 'post master' rank, based on input
176  // global index value
177  const std::int64_t num_nodes_g = mesh->geometry().index_map()->size_global();
178  const int comm_size = MPI::size(comm);
179  // NOTE: could make this int32_t be sending: index <- index - dest_rank_offset
180  std::vector<std::vector<std::int64_t>> nodes_g_send(comm_size);
181  for (std::int64_t node : nodes_g)
182  {
183  // TODO: Optimise this call by adding 'vectorised verion of
184  // MPI::index_owner
185  // Figure out which process is the postmaster for the input global index
186  const std::int32_t p
187  = dolfinx::MPI::index_owner(comm_size, node, num_nodes_g);
188  nodes_g_send[p].push_back(node);
189  }
190 
191  // Send/receive
192  const graph::AdjacencyList<std::int64_t> nodes_g_recv
194 
195  // -------------------
196  // 2. Send the entity key (nodes list) and tag to the postmaster based
197  // on the lowest index node in the entity 'key'
198  //
199  // NOTE: Stage 2 doesn't depend on the data received in Step 1, so
200  // data (i) the communication could be combined, or (ii) the
201  // communication in Step 1 could be make non-blocking.
202 
203  std::vector<std::vector<std::int64_t>> entities_send(comm_size);
204  std::vector<std::vector<T>> values_send(comm_size);
205  std::vector<std::int64_t> entity(entities.cols());
206  for (std::int32_t e = 0; e < entities.rows(); ++e)
207  {
208  // Copy nodes for entity and sort
209  std::copy(entities.row(e).data(), entities.row(e).data() + entities.cols(),
210  entity.begin());
211  std::sort(entity.begin(), entity.end());
212 
213  // Determine postmaster based on lowest entity node
214  const std::int32_t p
215  = dolfinx::MPI::index_owner(comm_size, entity.front(), num_nodes_g);
216  entities_send[p].insert(entities_send[p].end(), entity.begin(),
217  entity.end());
218  values_send[p].push_back(values[e]);
219  }
220 
221  // TODO: Pack into one MPI call
223  comm, graph::AdjacencyList<std::int64_t>(entities_send));
224  const graph::AdjacencyList<T> values_recv
225  = MPI::all_to_all(comm, graph::AdjacencyList<T>(values_send));
226 
227  // -------------------
228  // 3. As 'postmaster', send back the entity key (vertex list) and tag
229  // value to ranks that possibly need the data. Do this based on the
230  // first node index in the entity key.
231 
232  // NOTE: Could: (i) use a std::unordered_multimap, or (ii) only send
233  // owned nodes to the postmaster and use map, unordered_map or
234  // std::vector<pair>>, followed by a neighbourhood all_to_all at the
235  // end.
236  //
237  // Build map from global node index to ranks that have the node
238  std::multimap<std::int64_t, int> node_to_rank;
239  for (int p = 0; p < nodes_g_recv.num_nodes(); ++p)
240  {
241  auto nodes = nodes_g_recv.links(p);
242  for (int i = 0; i < nodes.rows(); ++i)
243  node_to_rank.insert({nodes(i), p});
244  }
245 
246  // Figure out which processes are owners of received nodes
247  std::vector<std::vector<std::int64_t>> send_nodes_owned(comm_size);
248  std::vector<std::vector<T>> send_vals_owned(comm_size);
249  const int nnodes_per_entity = entities.cols();
250  const Eigen::Map<const Eigen::Array<std::int64_t, Eigen::Dynamic,
251  Eigen::Dynamic, Eigen::RowMajor>>
252  _entities_recv(entities_recv.array().data(),
253  entities_recv.array().rows() / nnodes_per_entity,
254  nnodes_per_entity);
255  auto _values_recv = values_recv.array();
256  assert(_values_recv.rows() == _entities_recv.rows());
257  for (int e = 0; e < _entities_recv.rows(); ++e)
258  {
259  // Find ranks that have node0
260  auto [it0, it1] = node_to_rank.equal_range(_entities_recv(e, 0));
261  for (auto it = it0; it != it1; ++it)
262  {
263  const int p1 = it->second;
264  send_nodes_owned[p1].insert(
265  send_nodes_owned[p1].end(), _entities_recv.row(e).data(),
266  _entities_recv.row(e).data() + _entities_recv.cols());
267  send_vals_owned[p1].push_back(_values_recv(e));
268  }
269  }
270 
271  // TODO: Pack into one MPI call
273  comm, graph::AdjacencyList<std::int64_t>(send_nodes_owned));
274  const graph::AdjacencyList<T> recv_vals
275  = MPI::all_to_all(comm, graph::AdjacencyList<T>(send_vals_owned));
276 
277  // -------------------
278  // 4. From the received (key, value) data, determine which keys
279  // (entities) are on this process.
280 
281  // TODO: Rather than using std::map<std::vector<std::int64_t>,
282  // std::int32_t>, use a rectangular Eigen::Array to avoid the
283  // cost of std::vector<std::int64_t> allocations, and sort the
284  // Array by row.
285  //
286  // TODO: We have already received possibly tagged entities from other
287  // ranks, so we could use the received data to avoid creating
288  // the std::map for *all* entities and just for candidate
289  // entities.
290 
291  // Build map from vertex index (local to rank) to global "user" node
292  // index
293  auto map_v = mesh->topology().index_map(0);
294  assert(map_v);
295  const std::int32_t num_vertices = map_v->size_local() + map_v->num_ghosts();
296  const graph::AdjacencyList<std::int32_t>& x_dofmap
297  = mesh->geometry().dofmap();
298  std::vector<std::int64_t> vertex_to_node(num_vertices);
299  auto c_to_v = mesh->topology().connectivity(mesh->topology().dim(), 0);
300  if (!c_to_v)
301  throw std::runtime_error("missing cell-vertex connectivity.");
302  for (int c = 0; c < c_to_v->num_nodes(); ++c)
303  {
304  auto vertices = c_to_v->links(c);
305  auto x_dofs = x_dofmap.links(c);
306  for (int v = 0; v < vertices.rows(); ++v)
307  vertex_to_node[vertices[v]] = nodes_g[x_dofs[v]];
308  }
309 
310  // Build a map from entities on this process (keyed by vertex ordered
311  // entity input global indices) to entity local index
312  auto e_to_v = mesh->topology().connectivity(e_dim, 0);
313  if (!e_to_v)
314  throw std::runtime_error("Missing entity-vertex connectivity.");
315  std::map<std::vector<std::int64_t>, std::int32_t> entity_key_to_index;
316  std::vector<std::int64_t> key(nnodes_per_entity);
317  for (std::int32_t e = 0; e < e_to_v->num_nodes(); ++e)
318  {
319  auto vertices = e_to_v->links(e);
320  for (int v = 0; v < vertices.rows(); ++v)
321  key[v] = vertex_to_node[vertices(v)];
322  std::sort(key.begin(), key.end());
323  entity_key_to_index.insert({key, e});
324  }
325 
326  // Iterate over all received entities. If entity is on this rank,
327  // store (local entity index, tag value)
328  std::vector<std::int32_t> indices_new;
329  std::vector<T> values_new;
330  const Eigen::Map<const Eigen::Array<std::int64_t, Eigen::Dynamic,
331  Eigen::Dynamic, Eigen::RowMajor>>
332  _entities(recv_ents.array().data(),
333  recv_ents.array().rows() / nnodes_per_entity,
334  nnodes_per_entity);
335  entity.resize(nnodes_per_entity);
336  for (Eigen::Index e = 0; e < _entities.rows(); ++e)
337  {
338  // Note: _entities.row(e) was sorted by the sender
339  std::copy(_entities.row(e).data(),
340  _entities.row(e).data() + nnodes_per_entity, entity.begin());
341  if (const auto it = entity_key_to_index.find(entity);
342  it != entity_key_to_index.end())
343  {
344  indices_new.push_back(it->second);
345  values_new.push_back(recv_vals.array()[e]);
346  }
347  }
348 
349  // -------------------
350  // 5. Build MeshTags object
351 
352  auto [indices_sorted, values_sorted]
353  = common::sort_unique(indices_new, values_new);
354  return mesh::MeshTags<T>(mesh, e_dim, std::move(indices_sorted),
355  std::move(values_sorted));
356 }
357 } // namespace mesh
358 } // namespace dolfinx
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:109
dolfinx::mesh::CellType
CellType
Cell type identifier.
Definition: cell_types.h:22
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:103
dolfinx::mesh::create_meshtags
mesh::MeshTags< T > create_meshtags(MPI_Comm comm, const std::shared_ptr< const mesh::Mesh > &mesh, const mesh::CellType &tag_cell_type, const Eigen::Array< std::int64_t, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > &entities, const std::vector< T > &values)
Definition: MeshTags.h:147
dolfinx::mesh::cell_dim
int cell_dim(CellType type)
Return topological dimension of cell type.
Definition: cell_types.cpp:363
dolfinx::mesh::MeshTags::mesh
std::shared_ptr< const Mesh > mesh() const
Return mesh.
Definition: MeshTags.h:106
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::array
const Eigen::Array< T, Eigen::Dynamic, 1 > & array() const
Return contiguous array of links for all nodes (const version)
Definition: AdjacencyList.h:175
dolfinx::graph::AdjacencyList
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition: assemble_matrix_impl.h:26
dolfinx::MPI::size
static int size(MPI_Comm comm)
Return size of the group (number of processes) associated with the communicator.
Definition: MPI.cpp:83
dolfinx::mesh::MeshTags
A MeshTags are used to associate mesh entities with values. The entity index (local to process) ident...
Definition: Form.h:35
dolfinx::MPI::all_to_all
static graph::AdjacencyList< T > all_to_all(MPI_Comm comm, const graph::AdjacencyList< T > &send_data)
Send in_values[p0] to process p0 and receive values from process p1 in out_values[p1].
Definition: MPI.h:183
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:112
dolfinx::mesh::MeshTags::indices
const std::vector< std::int32_t > & indices() const
Indices of tagged mesh entities (local-to-process)
Definition: MeshTags.h:97
dolfinx::mesh::MeshTags::values
const std::vector< T > & values() const
Values attached to mesh entities.
Definition: MeshTags.h:100
dolfinx::MPI::index_owner
static int index_owner(int size, std::size_t index, std::size_t N)
Return which process owns index (inverse of local_range)
Definition: MPI.cpp:119
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.