OGS 6.2.0-97-g4a610c866
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  std::size_t offset = 0;
133  for (int variable_id = 0; variable_id < static_cast<int>(vec_var_n_components.size());
134  ++variable_id)
135  {
136  for (int component_id = 0; component_id < static_cast<int>(vec_var_n_components[variable_id]);
137  ++component_id)
138  {
139  auto const global_component_id =
140  getGlobalComponent(variable_id, component_id);
141 
142  auto const& ms = _mesh_subsets[global_component_id];
143  std::size_t const mesh_id = ms.getMeshID();
144 
145  findGlobalIndices(ms.elementsBegin(), ms.elementsEnd(),
146  ms.getNodes(), mesh_id, global_component_id,
147  global_component_id);
148  // increase by number of components of that variable
149  offset += _mesh_subsets.size();
150  }
151  }
152 }
153 
155  std::vector<MeshLib::MeshSubset>&& mesh_subsets,
156  std::vector<int> const& vec_var_n_components,
157  std::vector<std::vector<MeshLib::Element*> const*> const& vec_var_elements,
158  NumLib::ComponentOrder const order)
159  : _mesh_subsets(std::move(mesh_subsets)),
161  _variable_component_offsets(to_cumulative(vec_var_n_components))
162 {
163  assert(vec_var_n_components.size() == vec_var_elements.size());
164 
165  // For each element of that MeshSubset save a line of global indices.
166 
167  // _rows should be resized based on an element ID
168  std::size_t max_elem_id = 0;
169  for (std::vector<MeshLib::Element*>const* eles : vec_var_elements)
170  {
171  for (auto e : *eles)
172  {
173  max_elem_id = std::max(max_elem_id, e->getID());
174  }
175  }
176  _rows.resize(max_elem_id + 1, _mesh_subsets.size());
177 
178  std::size_t offset = 0;
179  for (int variable_id = 0; variable_id < static_cast<int>(vec_var_n_components.size());
180  ++variable_id)
181  {
182  std::vector<MeshLib::Element*> const& var_elements = *vec_var_elements[variable_id];
183  for (int component_id = 0; component_id < static_cast<int>(vec_var_n_components[variable_id]);
184  ++component_id)
185  {
186  auto const global_component_id =
187  getGlobalComponent(variable_id, component_id);
188 
189  auto const& ms = _mesh_subsets[global_component_id];
190  std::size_t const mesh_id = ms.getMeshID();
191 
193  var_elements.cbegin(), var_elements.cend(), ms.getNodes(),
194  mesh_id, global_component_id, global_component_id);
195  // increase by number of components of that variable
196  offset += _mesh_subsets.size();
197  }
198  }
199 }
200 
202  std::vector<MeshLib::MeshSubset>&& mesh_subsets,
203  std::vector<int> const& global_component_ids,
204  std::vector<int> const& variable_component_offsets,
205  std::vector<MeshLib::Element*> const& elements,
206  NumLib::MeshComponentMap&& mesh_component_map,
208  : _mesh_subsets(std::move(mesh_subsets)),
209  _mesh_component_map(std::move(mesh_component_map)),
210  _variable_component_offsets(variable_component_offsets)
211 {
212  // Each subset in the mesh_subsets represents a single component.
213  if (_mesh_subsets.size() != global_component_ids.size())
214  {
215  OGS_FATAL(
216  "Number of mesh subsets is not equal to number of components. "
217  "There are %d mesh subsets and %d components.",
218  _mesh_subsets.size(), global_component_ids.size());
219  }
220 
221  for (int i = 0; i < static_cast<int>(global_component_ids.size()); ++i)
222  {
223  auto const& ms = _mesh_subsets[i];
224 
225  // For all MeshSubset in mesh_subsets and each element of that
226  // MeshSubset save a line of global indices.
227  std::size_t const mesh_id = ms.getMeshID();
228 
229  findGlobalIndices(elements.cbegin(), elements.cend(), ms.getNodes(),
230  mesh_id, global_component_ids[i], i);
231  }
232 }
233 
235  int const variable_id,
236  std::vector<int> const& component_ids,
237  MeshLib::MeshSubset&& new_mesh_subset) const
238 {
239  DBUG("Construct reduced local to global index map.");
240 
241  if (component_ids.empty())
242  {
243  OGS_FATAL("Expected non-empty vector of component ids.");
244  }
245 
246  // Elements of the new_mesh_subset's mesh.
247  std::vector<MeshLib::Element*> const& elements =
248  new_mesh_subset.getMesh().getElements();
249 
250  // Create a subset of the current mesh component map.
251  std::vector<int> global_component_ids;
252 
253  for (auto component_id : component_ids)
254  {
255  global_component_ids.push_back(
256  getGlobalComponent(variable_id, component_id));
257  }
258 
259  auto mesh_component_map = _mesh_component_map.getSubset(
260  _mesh_subsets, new_mesh_subset, global_component_ids);
261 
262  // Create copies of the new_mesh_subset for each of the global components.
263  // The last component is moved after the for-loop.
264  std::vector<MeshLib::MeshSubset> all_mesh_subsets;
265  for (int i = 0; i < static_cast<int>(global_component_ids.size()) - 1; ++i)
266  {
267  all_mesh_subsets.emplace_back(new_mesh_subset);
268  }
269  all_mesh_subsets.emplace_back(std::move(new_mesh_subset));
270 
271  return new LocalToGlobalIndexMap(
272  std::move(all_mesh_subsets), global_component_ids,
273  _variable_component_offsets, elements, std::move(mesh_component_map),
274  ConstructorTag{});
275 }
276 
277 std::unique_ptr<LocalToGlobalIndexMap>
279  MeshLib::MeshSubset&& new_mesh_subset) const
280 {
281  DBUG("Construct reduced local to global index map.");
282 
283  // Create a subset of the current mesh component map.
284  std::vector<int> global_component_ids;
285 
286  for (int i = 0; i < getNumberOfComponents(); ++i)
287  {
288  global_component_ids.push_back(i);
289  }
290 
291  // Elements of the new_mesh_subset's mesh.
292  std::vector<MeshLib::Element*> const& elements =
293  new_mesh_subset.getMesh().getElements();
294 
295  auto mesh_component_map = _mesh_component_map.getSubset(
296  _mesh_subsets, new_mesh_subset, global_component_ids);
297 
298  // Create copies of the new_mesh_subset for each of the global components.
299  // The last component is moved after the for-loop.
300  std::vector<MeshLib::MeshSubset> all_mesh_subsets;
301  for (int i = 0; i < static_cast<int>(global_component_ids.size()) - 1; ++i)
302  {
303  all_mesh_subsets.emplace_back(new_mesh_subset);
304  }
305  all_mesh_subsets.emplace_back(std::move(new_mesh_subset));
306 
307  return std::make_unique<LocalToGlobalIndexMap>(
308  std::move(all_mesh_subsets), global_component_ids,
309  _variable_component_offsets, elements, std::move(mesh_component_map),
310  ConstructorTag{});
311 }
312 
314 {
316 }
317 
319 {
321 }
322 
324 {
325  return static_cast<int>(_variable_component_offsets.size()) - 1;
326 }
327 
329 {
330  assert(variable_id < getNumberOfVariables());
331  return _variable_component_offsets[variable_id + 1] -
332  _variable_component_offsets[variable_id];
333 }
334 
336 {
337  return _mesh_subsets.size();
338 }
339 
340 std::size_t LocalToGlobalIndexMap::size() const
341 {
342  return _rows.rows();
343 }
344 
346  std::size_t const mesh_item_id, const int component_id) const
347 {
348  return RowColumnIndices(_rows(mesh_item_id, component_id),
349  _columns(mesh_item_id, component_id));
350 }
351 
352 std::size_t
353 LocalToGlobalIndexMap::getNumberOfElementDOF(std::size_t const mesh_item_id) const
354 {
355  std::size_t ndof = 0;
356 
357  for (Table::Index c = 0; c < _rows.cols(); ++c)
358  {
359  ndof += _rows(mesh_item_id, c).size();
360  }
361 
362  return ndof;
363 }
364 
365 std::size_t
366 LocalToGlobalIndexMap::getNumberOfElementComponents(std::size_t const mesh_item_id) const
367 {
368  std::size_t n = 0;
369  for (Table::Index c = 0; c < _rows.cols(); ++c)
370  {
371  if (!_rows(mesh_item_id, c).empty())
372  {
373  n++;
374  }
375  }
376  return n;
377 }
378 
380  std::size_t const mesh_item_id) const
381 {
382  std::vector<int> vec;
383  for (int i = 0; i < getNumberOfVariables(); i++)
384  {
385  for (int j=0; j<getNumberOfVariableComponents(i); j++)
386  {
387  auto comp_id = getGlobalComponent(i, j);
388  if (!_rows(mesh_item_id, comp_id).empty())
389  {
390  vec.push_back(i);
391  }
392  }
393  }
394  std::sort(vec.begin(), vec.end());
395  vec.erase(std::unique(vec.begin(), vec.end()), vec.end());
396 
397  return vec;
398 }
399 
401  MeshLib::Location const& l,
402  int const variable_id,
403  int const component_id) const
404 {
405  auto const c = getGlobalComponent(variable_id, component_id);
407 }
408 
410  MeshLib::Location const& l, int const global_component_id) const
411 {
412  return _mesh_component_map.getGlobalIndex(l, global_component_id);
413 }
414 
416 std::vector<GlobalIndexType> LocalToGlobalIndexMap::getGlobalIndices(
417  const MeshLib::Location& l) const
418 {
420 }
421 
423 std::vector<GlobalIndexType> const& LocalToGlobalIndexMap::getGhostIndices()
424  const
425 {
427 }
428 
432  MeshLib::Location const& l, std::size_t const comp_id,
433  std::size_t const range_begin, std::size_t const range_end) const
434 {
435  return _mesh_component_map.getLocalIndex(l, comp_id, range_begin,
436  range_end);
437 }
438 
440  int const variable_id, int const component_id) const
441 {
442  return getMeshSubset(getGlobalComponent(variable_id, component_id));
443 }
444 
446  int const global_component_id) const
447 {
448  return _mesh_subsets[global_component_id];
449 }
450 
451 #ifndef NDEBUG
452 std::ostream& operator<<(std::ostream& os, LocalToGlobalIndexMap const& map)
453 {
454  std::size_t const max_lines = 10;
455  std::size_t lines_printed = 0;
456 
457  os << "Rows of the local to global index map; " << map._rows.size()
458  << " rows\n";
459  for (std::size_t e=0; e<map.size(); ++e)
460  {
461  os << "== e " << e << " ==\n";
462  for (int c = 0; c < map.getNumberOfComponents(); ++c)
463  {
464  auto const& line = map._rows(e, c);
465 
466  os << "c" << c << " { ";
467  std::copy(line.cbegin(), line.cend(),
468  std::ostream_iterator<std::size_t>(os, " "));
469  os << " }\n";
470  }
471 
472  if (lines_printed++ > max_lines)
473  {
474  os << "...\n";
475  break;
476  }
477  }
478  lines_printed = 0;
479 
480  os << "Mesh component map:\n" << map._mesh_component_map;
481  return os;
482 }
483 #endif // NDEBUG
484 
485 } // 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
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
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
void findGlobalIndices(ElementIterator first, ElementIterator last, std::vector< MeshLib::Node *> const &nodes, std::size_t const mesh_id, const int component_id, const int comp_id_write)
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.
void findGlobalIndicesWithElementID(ElementIterator first, ElementIterator last, std::vector< MeshLib::Node *> const &nodes, std::size_t const mesh_id, const int component_id, const int comp_id_write)
#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