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