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>
46 template <
typename U,
typename V>
50 _values(std::forward<V>(
values))
54 throw std::runtime_error(
55 "Indices and values arrays must have same size.");
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");
83 const Eigen::Array<std::int32_t, Eigen::Dynamic, 1>
find(
const T value)
const
85 int n = std::count(_values.begin(), _values.end(), value);
86 Eigen::Array<std::int32_t, Eigen::Dynamic, 1>
indices(n);
88 for (std::int32_t i = 0; i < _values.size(); ++i)
90 if (_values[i] == value)
91 indices[counter++] = _indices[i];
97 const std::vector<std::int32_t>&
indices()
const {
return _indices; }
100 const std::vector<T>&
values()
const {
return _values; }
103 int dim()
const {
return _dim; }
106 std::shared_ptr<const Mesh>
mesh()
const {
return _mesh; }
109 std::string
name =
"mesh_tags";
112 std::size_t
id()
const {
return _unique_id; }
119 std::shared_ptr<const Mesh> _mesh;
125 std::vector<std::int32_t> _indices;
128 std::vector<T> _values;
145 template <
typename T>
149 const Eigen::Array<std::int64_t, Eigen::Dynamic, Eigen::Dynamic,
150 Eigen::RowMajor>& entities,
151 const std::vector<T>& values)
159 if ((std::size_t)entities.rows() != values.size())
160 throw std::runtime_error(
"Number of entities and values must match");
172 const std::vector<std::int64_t>& nodes_g
173 = mesh->geometry().input_global_indices();
177 const std::int64_t num_nodes_g = mesh->geometry().index_map()->size_global();
180 std::vector<std::vector<std::int64_t>> nodes_g_send(comm_size);
181 for (std::int64_t node : nodes_g)
188 nodes_g_send[p].push_back(node);
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)
209 std::copy(entities.row(e).data(), entities.row(e).data() + entities.cols(),
211 std::sort(entity.begin(), entity.end());
216 entities_send[p].insert(entities_send[p].end(), entity.begin(),
218 values_send[p].push_back(values[e]);
238 std::multimap<std::int64_t, int> node_to_rank;
239 for (
int p = 0; p < nodes_g_recv.
num_nodes(); ++p)
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});
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,
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)
260 auto [it0, it1] = node_to_rank.equal_range(_entities_recv(e, 0));
261 for (
auto it = it0; it != it1; ++it)
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));
293 auto map_v = mesh->topology().index_map(0);
295 const std::int32_t num_vertices = map_v->size_local() + map_v->num_ghosts();
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);
301 throw std::runtime_error(
"missing cell-vertex connectivity.");
302 for (
int c = 0; c < c_to_v->num_nodes(); ++c)
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]];
312 auto e_to_v = mesh->topology().connectivity(e_dim, 0);
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)
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});
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,
335 entity.resize(nnodes_per_entity);
336 for (Eigen::Index e = 0; e < _entities.rows(); ++e)
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())
344 indices_new.push_back(it->second);
345 values_new.push_back(recv_vals.
array()[e]);
352 auto [indices_sorted, values_sorted]
355 std::move(values_sorted));