OGS
LocalToGlobalIndexMap.cpp
Go to the documentation of this file.
1
12
13#include <algorithm>
14#include <numeric>
15#include <unordered_set>
16
17namespace NumLib
18{
19namespace
20{
21// Make the cumulative sum of an array, which starts with zero
22template <typename T>
23std::vector<T> to_cumulative(std::vector<T> const& vec)
24{
25 std::vector<T> result(vec.size() + 1, 0);
26 std::partial_sum(vec.begin(), vec.end(), result.begin() + 1);
27
28 return result;
29}
30
31} // namespace
32
34 int const component_id) const
35{
36 return _variable_component_offsets[variable_id] + component_id;
37}
38
39template <typename ElementIterator>
41 ElementIterator first, ElementIterator last,
42 std::vector<MeshLib::Node*> const& nodes, std::size_t const mesh_id,
43 const int comp_id, const int comp_id_write)
44{
45 std::unordered_set<MeshLib::Node*> const set_nodes(nodes.begin(),
46 nodes.end());
47
48 // For each element find the global indices for node/element
49 // components.
50 for (ElementIterator e = first; e != last; ++e)
51 {
52 LineIndex indices;
53 indices.reserve((*e)->getNumberOfNodes());
54
55 for (auto* n = (*e)->getNodes();
56 n < (*e)->getNodes() + (*e)->getNumberOfNodes();
57 ++n)
58 {
59 // Check if the element's node is in the given list of nodes.
60 if (set_nodes.find(*n) == set_nodes.end())
61 {
62 continue;
63 }
65 (*n)->getID());
66 indices.push_back(_mesh_component_map.getGlobalIndex(l, comp_id));
67 }
68
69 indices.shrink_to_fit();
70 _rows((*e)->getID(), comp_id_write) = std::move(indices);
71 }
72}
73
74template <typename ElementIterator>
76 ElementIterator first, ElementIterator last,
77 std::vector<MeshLib::Node*> const& nodes, std::size_t const mesh_id,
78 const int comp_id, const int comp_id_write)
79{
80 _rows.resize(std::distance(first, last), _mesh_subsets.size());
81
82 std::unordered_set<MeshLib::Node*> const set_nodes(nodes.begin(),
83 nodes.end());
84
85 // For each element find the global indices for node/element
86 // components.
87 std::size_t elem_id = 0;
88 for (ElementIterator e = first; e != last; ++e, ++elem_id)
89 {
90 LineIndex indices;
91 indices.reserve((*e)->getNumberOfNodes());
92
93 for (auto* n = (*e)->getNodes();
94 n < (*e)->getNodes() + (*e)->getNumberOfNodes();
95 ++n)
96 {
97 // Check if the element's node is in the given list of nodes.
98 if (set_nodes.find(*n) == set_nodes.end())
99 {
100 continue;
101 }
103 (*n)->getID());
104 auto const global_index =
106 if (global_index == std::numeric_limits<GlobalIndexType>::max())
107 {
108 continue;
109 }
110 indices.push_back(global_index);
111 }
112
113 indices.shrink_to_fit();
114 _rows(elem_id, comp_id_write) = std::move(indices);
115 }
116}
117
119 std::vector<MeshLib::MeshSubset>&& mesh_subsets,
120 NumLib::ComponentOrder const order)
121 : LocalToGlobalIndexMap(std::move(mesh_subsets),
122 std::vector<int>(mesh_subsets.size(), 1), order)
123{
124}
125
127 std::vector<MeshLib::MeshSubset>&& mesh_subsets,
128 std::vector<int> const& vec_var_n_components,
129 NumLib::ComponentOrder const order)
130 : _mesh_subsets(std::move(mesh_subsets)),
131 _mesh_component_map(_mesh_subsets, order),
132 _variable_component_offsets(to_cumulative(vec_var_n_components))
133{
134 // For each element of that MeshSubset save a line of global indices.
135 for (int variable_id = 0;
136 variable_id < static_cast<int>(vec_var_n_components.size());
137 ++variable_id)
138 {
139 for (int component_id = 0;
140 component_id < static_cast<int>(vec_var_n_components[variable_id]);
141 ++component_id)
142 {
143 auto const global_component_id =
144 getGlobalComponent(variable_id, component_id);
145
146 auto const& ms = _mesh_subsets[global_component_id];
147 std::size_t const mesh_id = ms.getMeshID();
148
149 findGlobalIndices(ms.elementsBegin(), ms.elementsEnd(),
150 ms.getNodes(), mesh_id, global_component_id,
151 global_component_id);
152 }
153 }
154}
155
157 std::vector<MeshLib::MeshSubset>&& mesh_subsets,
158 std::vector<int> const& vec_var_n_components,
159 std::vector<std::vector<MeshLib::Element*> const*> const& vec_var_elements,
160 NumLib::ComponentOrder const order)
161 : _mesh_subsets(std::move(mesh_subsets)),
162 _mesh_component_map(_mesh_subsets, order),
163 _variable_component_offsets(to_cumulative(vec_var_n_components))
164{
165 assert(vec_var_n_components.size() == vec_var_elements.size());
166
167 // For each element of that MeshSubset save a line of global indices.
168
169 // _rows should be resized based on an element ID
170 std::size_t max_elem_id = 0;
171 for (std::vector<MeshLib::Element*> const* elements : vec_var_elements)
172 {
173 for (auto e : *elements)
174 {
175 max_elem_id = std::max(max_elem_id, e->getID());
176 }
177 }
178 _rows.resize(max_elem_id + 1, _mesh_subsets.size());
179
180 for (int variable_id = 0;
181 variable_id < static_cast<int>(vec_var_n_components.size());
182 ++variable_id)
183 {
184 std::vector<MeshLib::Element*> const& var_elements =
185 *vec_var_elements[variable_id];
186 for (int component_id = 0;
187 component_id < static_cast<int>(vec_var_n_components[variable_id]);
188 ++component_id)
189 {
190 auto const global_component_id =
191 getGlobalComponent(variable_id, component_id);
192
193 auto const& ms = _mesh_subsets[global_component_id];
194 std::size_t const mesh_id = ms.getMeshID();
195
197 var_elements.cbegin(), var_elements.cend(), ms.getNodes(),
198 mesh_id, global_component_id, global_component_id);
199 }
200 }
201}
202
204 std::vector<MeshLib::MeshSubset>&& mesh_subsets,
205 std::vector<int> const& global_component_ids,
206 std::vector<int> const& variable_component_offsets,
207 std::vector<MeshLib::Element*> const& elements,
208 NumLib::MeshComponentMap&& mesh_component_map,
210 : _mesh_subsets(std::move(mesh_subsets)),
211 _mesh_component_map(std::move(mesh_component_map)),
212 _variable_component_offsets(variable_component_offsets)
213{
214 // Each subset in the mesh_subsets represents a single component.
215 if (_mesh_subsets.size() != global_component_ids.size())
216 {
217 OGS_FATAL(
218 "Number of mesh subsets is not equal to number of components. "
219 "There are {:d} mesh subsets and {:d} components.",
220 _mesh_subsets.size(), global_component_ids.size());
221 }
222
223 for (int i = 0; i < static_cast<int>(global_component_ids.size()); ++i)
224 {
225 auto const& ms = _mesh_subsets[i];
226
227 // For all MeshSubset in mesh_subsets and each element of that
228 // MeshSubset save a line of global indices.
229 std::size_t const mesh_id = ms.getMeshID();
230
231 findGlobalIndices(elements.cbegin(), elements.cend(), ms.getNodes(),
232 mesh_id, global_component_ids[i], i);
233 }
234}
235
236std::unique_ptr<LocalToGlobalIndexMap>
238 int const variable_id,
239 std::vector<int> const& component_ids,
240 MeshLib::MeshSubset&& new_mesh_subset) const
241{
242 DBUG("Construct reduced local to global index map.");
243
244 if (component_ids.empty())
245 {
246 OGS_FATAL("Expected non-empty vector of component ids.");
247 }
248
249 // Elements of the new_mesh_subset's mesh.
250 std::vector<MeshLib::Element*> const& elements =
251 new_mesh_subset.getMesh().getElements();
252
253 // Create a subset of the current mesh component map.
254 std::vector<int> global_component_ids;
255
256 transform(cbegin(component_ids), cend(component_ids),
257 back_inserter(global_component_ids),
258 [&](auto const component_id)
259 { return getGlobalComponent(variable_id, component_id); });
260
261 auto mesh_component_map = _mesh_component_map.getSubset(
262 _mesh_subsets, new_mesh_subset, global_component_ids);
263
264 // Create copies of the new_mesh_subset for each of the global components.
265 // The last component is moved after the for-loop.
266 std::vector<MeshLib::MeshSubset> all_mesh_subsets;
267 for (int i = 0; i < static_cast<int>(global_component_ids.size()) - 1; ++i)
268 {
269 all_mesh_subsets.emplace_back(new_mesh_subset);
270 }
271 all_mesh_subsets.emplace_back(std::move(new_mesh_subset));
272
273 return std::make_unique<LocalToGlobalIndexMap>(
274 std::move(all_mesh_subsets), global_component_ids,
275 _variable_component_offsets, elements, std::move(mesh_component_map),
277}
278
279std::unique_ptr<LocalToGlobalIndexMap>
281 MeshLib::MeshSubset&& new_mesh_subset) const
282{
283 DBUG("Construct reduced local to global index map.");
284
285 // Create a subset of the current mesh component map.
286 std::vector<int> global_component_ids;
287
288 for (int i = 0; i < getNumberOfGlobalComponents(); ++i)
289 {
290 global_component_ids.push_back(i);
291 }
292
293 // Elements of the new_mesh_subset's mesh.
294 std::vector<MeshLib::Element*> const& elements =
295 new_mesh_subset.getMesh().getElements();
296
297 auto mesh_component_map = _mesh_component_map.getSubset(
298 _mesh_subsets, new_mesh_subset, global_component_ids);
299
300 // Create copies of the new_mesh_subset for each of the global components.
301 // The last component is moved after the for-loop.
302 std::vector<MeshLib::MeshSubset> all_mesh_subsets;
303 for (int i = 0; i < static_cast<int>(global_component_ids.size()) - 1; ++i)
304 {
305 all_mesh_subsets.emplace_back(new_mesh_subset);
306 }
307 all_mesh_subsets.emplace_back(std::move(new_mesh_subset));
308
309 return std::make_unique<LocalToGlobalIndexMap>(
310 std::move(all_mesh_subsets), global_component_ids,
311 _variable_component_offsets, elements, std::move(mesh_component_map),
313}
314
319
324
326{
327 return static_cast<int>(_variable_component_offsets.size()) - 1;
328}
329
331{
332 assert(variable_id < getNumberOfVariables());
333 return _variable_component_offsets[variable_id + 1] -
334 _variable_component_offsets[variable_id];
335}
336
341
343{
344 return _rows.rows();
345}
346
348 std::size_t const mesh_item_id, const int global_component_id) const
349{
350 return RowColumnIndices(_rows(mesh_item_id, global_component_id),
351 _columns(mesh_item_id, global_component_id));
352}
353
355 std::size_t const mesh_item_id) const
356{
357 std::size_t ndof = 0;
358
359 for (Table::Index c = 0; c < _rows.cols(); ++c)
360 {
361 ndof += _rows(mesh_item_id, c).size();
362 }
363
364 return ndof;
365}
366
368 std::size_t const mesh_item_id) const
369{
370 std::size_t n = 0;
371 for (Table::Index c = 0; c < _rows.cols(); ++c)
372 {
373 if (!_rows(mesh_item_id, c).empty())
374 {
375 n++;
376 }
377 }
378 return n;
379}
380
382 std::size_t const mesh_item_id) const
383{
384 std::vector<int> vec;
385 for (int i = 0; i < getNumberOfVariables(); i++)
386 {
387 for (int j = 0; j < getNumberOfVariableComponents(i); j++)
388 {
389 auto comp_id = getGlobalComponent(i, j);
390 if (!_rows(mesh_item_id, comp_id).empty())
391 {
392 vec.push_back(i);
393 }
394 }
395 }
396 std::sort(vec.begin(), vec.end());
397 vec.erase(std::unique(vec.begin(), vec.end()), vec.end());
398
399 return vec;
400}
401
403 MeshLib::Location const& l,
404 int const variable_id,
405 int const component_id) const
406{
407 auto const c = getGlobalComponent(variable_id, component_id);
409}
410
412 MeshLib::Location const& l, int const global_component_id) const
413{
414 return _mesh_component_map.getGlobalIndex(l, global_component_id);
415}
416
418std::vector<GlobalIndexType> LocalToGlobalIndexMap::getGlobalIndices(
419 const MeshLib::Location& l) const
420{
422}
423
425std::vector<GlobalIndexType> const& LocalToGlobalIndexMap::getGhostIndices()
426 const
427{
429}
430
434 MeshLib::Location const& l, std::size_t const comp_id,
435 std::size_t const range_begin, std::size_t const range_end) const
436{
437 return _mesh_component_map.getLocalIndex(l, comp_id, range_begin,
438 range_end);
439}
440
442 int const variable_id, int const component_id) const
443{
444 return getMeshSubset(getGlobalComponent(variable_id, component_id));
445}
446
448 int const global_component_id) const
449{
450 return _mesh_subsets[global_component_id];
451}
452
453#ifndef NDEBUG
454std::ostream& operator<<(std::ostream& os, LocalToGlobalIndexMap const& map)
455{
456 std::size_t const max_lines = 10;
457 std::size_t lines_printed = 0;
458
459 os << "Rows of the local to global index map; " << map._rows.size()
460 << " rows\n";
461 for (std::size_t e = 0; e < map.size(); ++e)
462 {
463 os << "== e " << e << " ==\n";
464 for (int c = 0; c < map.getNumberOfGlobalComponents(); ++c)
465 {
466 auto const& line = map._rows(e, c);
467
468 os << "c" << c << " { ";
469 std::copy(line.cbegin(), line.cend(),
470 std::ostream_iterator<std::size_t>(os, " "));
471 os << " }\n";
472 }
473
474 if (lines_printed++ > max_lines)
475 {
476 os << "...\n";
477 break;
478 }
479 }
480
481 os << "Mesh component map:\n" << map._mesh_component_map;
482 return os;
483}
484#endif // NDEBUG
485
486} // namespace NumLib
#define OGS_FATAL(...)
Definition Error.h:26
GlobalMatrix::IndexType GlobalIndexType
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
A subset of nodes on a single mesh.
Definition MeshSubset.h:26
RowColumnIndices operator()(std::size_t const mesh_item_id, const int global_component_id) const
std::size_t getNumberOfElementDOF(std::size_t const mesh_item_id) const
LocalToGlobalIndexMap(std::vector< MeshLib::MeshSubset > &&mesh_subsets, NumLib::ComponentOrder const order)
std::vector< MeshLib::MeshSubset > _mesh_subsets
A vector of mesh subsets for each process variables' components.
int getGlobalComponent(int const variable_id, int const component_id) const
NumLib::MeshComponentMap _mesh_component_map
MathLib::RowColumnIndices< GlobalIndexType > RowColumnIndices
GlobalIndexType getGlobalIndex(MeshLib::Location const &l, int const variable_id, int const component_id) const
void findGlobalIndicesWithElementID(ElementIterator first, ElementIterator last, std::vector< MeshLib::Node * > const &nodes, std::size_t const mesh_id, const int comp_id, const int comp_id_write)
std::vector< int > const _variable_component_offsets
int getNumberOfVariableComponents(int variable_id) const
std::unique_ptr< LocalToGlobalIndexMap > deriveBoundaryConstrainedMap(int const variable_id, std::vector< int > const &component_ids, MeshLib::MeshSubset &&new_mesh_subset) const
MeshLib::MeshSubset const & getMeshSubset(int const variable_id, int const component_id) const
std::size_t getNumberOfElementComponents(std::size_t const mesh_item_id) const
void findGlobalIndices(ElementIterator first, ElementIterator last, std::vector< MeshLib::Node * > const &nodes, std::size_t const mesh_id, const int comp_id, const int comp_id_write)
GlobalIndexType getLocalIndex(MeshLib::Location const &l, std::size_t const comp_id, std::size_t const range_begin, std::size_t const range_end) const
std::vector< GlobalIndexType > const & getGhostIndices() const
Get ghost indices, forwarded from MeshComponentMap.
RowColumnIndices::LineIndex LineIndex
std::vector< GlobalIndexType > getGlobalIndices(const MeshLib::Location &l) const
Forwards the respective method from MeshComponentMap.
std::vector< int > getElementVariableIDs(std::size_t const mesh_item_id) const
Multidirectional mapping between mesh entities and degrees of freedom.
std::size_t dofSizeWithoutGhosts() const
GlobalIndexType getGlobalIndex(Location const &l, int const comp_id) const
GlobalIndexType getLocalIndex(Location const &l, int const comp_id, std::size_t const range_begin, std::size_t const range_end) const
std::vector< GlobalIndexType > getGlobalIndices(const Location &l) const
MeshComponentMap getSubset(std::vector< MeshLib::MeshSubset > const &bulk_mesh_subsets, MeshLib::MeshSubset const &new_mesh_subset, std::vector< int > const &new_global_component_ids) const
std::size_t dofSizeWithGhosts() const
The number of dofs including the those located in the ghost nodes.
std::vector< GlobalIndexType > const & getGhostIndices() const
Get ghost indices (for DDC).
std::vector< T > to_cumulative(std::vector< T > const &vec)
ComponentOrder
Ordering of components in global matrix/vector.
std::ostream & operator<<(std::ostream &os, LocalToGlobalIndexMap const &map)