OGS
LocalToGlobalIndexMap.cpp
Go to the documentation of this file.
1 
11 #include "LocalToGlobalIndexMap.h"
12 
13 #include <algorithm>
14 #include <numeric>
15 #include <unordered_set>
16 
17 namespace NumLib
18 {
19 namespace
20 {
21 // Make the cumulative sum of an array, which starts with zero
22 template <typename T>
23 std::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 
39 template <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 
74 template <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 
236 std::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),
276  ConstructorTag{});
277 }
278 
279 std::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),
312  ConstructorTag{});
313 }
314 
316 {
318 }
319 
321 {
323 }
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 
338 {
339  return _mesh_subsets.size();
340 }
341 
342 std::size_t LocalToGlobalIndexMap::size() const
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 
418 std::vector<GlobalIndexType> LocalToGlobalIndexMap::getGlobalIndices(
419  const MeshLib::Location& l) const
420 {
422 }
423 
425 std::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
454 std::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(char const *fmt, Args const &... args)
Definition: Logging.h:27
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
std::vector< GlobalIndexType > const & getGhostIndices() const
Get ghost indices (for DDC).
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.
void copy(PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:37
std::vector< T > to_cumulative(std::vector< T > const &vec)
std::ostream & operator<<(std::ostream &os, LocalToGlobalIndexMap const &map)
ComponentOrder
Ordering of components in global matrix/vector.