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];
98 const std::vector<std::int32_t>&
indices()
const {
return _indices; }
101 const std::vector<T>&
values()
const {
return _values; }
104 int dim()
const {
return _dim; }
107 std::shared_ptr<const Mesh>
mesh()
const {
return _mesh; }
110 std::string
name =
"mesh_tags";
113 std::size_t
id()
const {
return _unique_id; }
120 std::shared_ptr<const Mesh> _mesh;
126 std::vector<std::int32_t> _indices;
129 std::vector<T> _values;
138 template <
typename T>
142 const std::vector<T>& values)
145 if ((std::size_t)entities.
num_nodes() != values.size())
146 throw std::runtime_error(
"Number of entities and values must match");
149 const auto map_e = mesh->topology().index_map(dim);
152 auto e_to_v = mesh->topology().connectivity(dim, 0);
154 throw std::runtime_error(
"Missing entity-vertex connectivity.");
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();
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)
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});
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);
178 for (Eigen::Index e = 0; e < entities.
num_nodes(); ++e)
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,
185 std::sort(entity.begin(), entity.end());
187 if (
const auto it = entity_key_to_index.find(entity);
188 it != entity_key_to_index.end())
190 indices_new.push_back(it->second);
191 values_new.push_back(values[e]);
195 auto [indices_sorted, values_sorted]
198 std::move(values_sorted));