OGS 6.1.0-1721-g6382411ad
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 } // no named 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  continue;
62  mesh_id, MeshLib::MeshItemType::Node, (*n)->getID());
63  indices.push_back(_mesh_component_map.getGlobalIndex(l, comp_id));
64  }
65 
66  indices.shrink_to_fit();
67  _rows((*e)->getID(), comp_id_write) = std::move(indices);
68  }
69 }
70 
71 template <typename ElementIterator>
73  ElementIterator first, ElementIterator last,
74  std::vector<MeshLib::Node*> const& nodes, std::size_t const mesh_id,
75  const int comp_id, const int comp_id_write)
76 {
77  _rows.resize(std::distance(first, last), _mesh_subsets.size());
78 
79  std::unordered_set<MeshLib::Node*> const set_nodes(nodes.begin(), nodes.end());
80 
81  // For each element find the global indices for node/element
82  // components.
83  std::size_t elem_id = 0;
84  for (ElementIterator e = first; e != last; ++e, ++elem_id)
85  {
86  LineIndex indices;
87  indices.reserve((*e)->getNumberOfNodes());
88 
89  for (auto* n = (*e)->getNodes();
90  n < (*e)->getNodes() + (*e)->getNumberOfNodes(); ++n)
91  {
92  // Check if the element's node is in the given list of nodes.
93  if (set_nodes.find(*n)==set_nodes.end())
94  continue;
96  mesh_id, MeshLib::MeshItemType::Node, (*n)->getID());
97  auto const global_index =
99  if (global_index == std::numeric_limits<GlobalIndexType>::max())
100  {
101  continue;
102  }
103  indices.push_back(global_index);
104  }
105 
106  indices.shrink_to_fit();
107  _rows(elem_id, comp_id_write) = std::move(indices);
108  }
109 }
110 
112  std::vector<MeshLib::MeshSubset>&& mesh_subsets,
113  NumLib::ComponentOrder const order)
114  : LocalToGlobalIndexMap(std::move(mesh_subsets),
115  std::vector<int>(mesh_subsets.size(), 1), order)
116 {
117 }
118 
120  std::vector<MeshLib::MeshSubset>&& mesh_subsets,
121  std::vector<int> const& vec_var_n_components,
122  NumLib::ComponentOrder const order)
123  : _mesh_subsets(std::move(mesh_subsets)),
125  _variable_component_offsets(to_cumulative(vec_var_n_components))
126 {
127  // For each element of that MeshSubset save a line of global indices.
128  std::size_t offset = 0;
129  for (int variable_id = 0; variable_id < static_cast<int>(vec_var_n_components.size());
130  ++variable_id)
131  {
132  for (int component_id = 0; component_id < static_cast<int>(vec_var_n_components[variable_id]);
133  ++component_id)
134  {
135  auto const global_component_id =
136  getGlobalComponent(variable_id, component_id);
137 
138  auto const& ms = _mesh_subsets[global_component_id];
139  std::size_t const mesh_id = ms.getMeshID();
140 
141  findGlobalIndices(ms.elementsBegin(), ms.elementsEnd(),
142  ms.getNodes(), mesh_id, global_component_id,
143  global_component_id);
144  // increase by number of components of that variable
145  offset += _mesh_subsets.size();
146  }
147  }
148 }
149 
151  std::vector<MeshLib::MeshSubset>&& mesh_subsets,
152  std::vector<int> const& vec_var_n_components,
153  std::vector<std::vector<MeshLib::Element*> const*> const& vec_var_elements,
154  NumLib::ComponentOrder const order)
155  : _mesh_subsets(std::move(mesh_subsets)),
157  _variable_component_offsets(to_cumulative(vec_var_n_components))
158 {
159  assert(vec_var_n_components.size() == vec_var_elements.size());
160 
161  // For each element of that MeshSubset save a line of global indices.
162 
163  // _rows should be resized based on an element ID
164  std::size_t max_elem_id = 0;
165  for (std::vector<MeshLib::Element*>const* eles : vec_var_elements)
166  {
167  for (auto e : *eles)
168  max_elem_id = std::max(max_elem_id, e->getID());
169  }
170  _rows.resize(max_elem_id + 1, _mesh_subsets.size());
171 
172  std::size_t offset = 0;
173  for (int variable_id = 0; variable_id < static_cast<int>(vec_var_n_components.size());
174  ++variable_id)
175  {
176  std::vector<MeshLib::Element*> const& var_elements = *vec_var_elements[variable_id];
177  for (int component_id = 0; component_id < static_cast<int>(vec_var_n_components[variable_id]);
178  ++component_id)
179  {
180  auto const global_component_id =
181  getGlobalComponent(variable_id, component_id);
182 
183  auto const& ms = _mesh_subsets[global_component_id];
184  std::size_t const mesh_id = ms.getMeshID();
185 
187  var_elements.cbegin(), var_elements.cend(), ms.getNodes(),
188  mesh_id, global_component_id, global_component_id);
189  // increase by number of components of that variable
190  offset += _mesh_subsets.size();
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  OGS_FATAL(
209  "Number of mesh subsets is not equal to number of components. "
210  "There are %d mesh subsets and %d components.",
211  _mesh_subsets.size(), global_component_ids.size());
212 
213  for (int i = 0; i < static_cast<int>(global_component_ids.size()); ++i)
214  {
215  auto const& ms = _mesh_subsets[i];
216 
217  // For all MeshSubset in mesh_subsets and each element of that
218  // MeshSubset save a line of global indices.
219  std::size_t const mesh_id = ms.getMeshID();
220 
221  findGlobalIndices(elements.cbegin(), elements.cend(), ms.getNodes(),
222  mesh_id, global_component_ids[i], i);
223  }
224 }
225 
227  int const variable_id,
228  std::vector<int> const& component_ids,
229  MeshLib::MeshSubset&& new_mesh_subset) const
230 {
231  DBUG("Construct reduced local to global index map.");
232 
233  if (component_ids.empty())
234  OGS_FATAL("Expected non-empty vector of component ids.");
235 
236  // Elements of the new_mesh_subset's mesh.
237  std::vector<MeshLib::Element*> const& elements =
238  new_mesh_subset.getMesh().getElements();
239 
240  // Create a subset of the current mesh component map.
241  std::vector<int> global_component_ids;
242 
243  for (auto component_id : component_ids)
244  global_component_ids.push_back(
245  getGlobalComponent(variable_id, component_id));
246 
247  auto mesh_component_map = _mesh_component_map.getSubset(
248  _mesh_subsets, new_mesh_subset, global_component_ids);
249 
250  // Create copies of the new_mesh_subset for each of the global components.
251  // The last component is moved after the for-loop.
252  std::vector<MeshLib::MeshSubset> all_mesh_subsets;
253  for (int i = 0; i < static_cast<int>(global_component_ids.size()) - 1; ++i)
254  all_mesh_subsets.emplace_back(new_mesh_subset);
255  all_mesh_subsets.emplace_back(std::move(new_mesh_subset));
256 
257  return new LocalToGlobalIndexMap(
258  std::move(all_mesh_subsets), global_component_ids,
259  _variable_component_offsets, elements, std::move(mesh_component_map),
260  ConstructorTag{});
261 }
262 
263 std::unique_ptr<LocalToGlobalIndexMap>
265  MeshLib::MeshSubset&& new_mesh_subset) const
266 {
267  DBUG("Construct reduced local to global index map.");
268 
269  // Create a subset of the current mesh component map.
270  std::vector<int> global_component_ids;
271 
272  for (int i = 0; i < getNumberOfComponents(); ++i)
273  {
274  global_component_ids.push_back(i);
275  }
276 
277  // Elements of the new_mesh_subset's mesh.
278  std::vector<MeshLib::Element*> const& elements =
279  new_mesh_subset.getMesh().getElements();
280 
281  auto mesh_component_map = _mesh_component_map.getSubset(
282  _mesh_subsets, new_mesh_subset, global_component_ids);
283 
284  // Create copies of the new_mesh_subset for each of the global components.
285  // The last component is moved after the for-loop.
286  std::vector<MeshLib::MeshSubset> all_mesh_subsets;
287  for (int i = 0; i < static_cast<int>(global_component_ids.size()) - 1; ++i)
288  all_mesh_subsets.emplace_back(new_mesh_subset);
289  all_mesh_subsets.emplace_back(std::move(new_mesh_subset));
290 
291  return std::make_unique<LocalToGlobalIndexMap>(
292  std::move(all_mesh_subsets), global_component_ids,
293  _variable_component_offsets, elements, std::move(mesh_component_map),
294  ConstructorTag{});
295 }
296 
298 {
300 }
301 
303 {
305 }
306 
308 {
309  return static_cast<int>(_variable_component_offsets.size()) - 1;
310 }
311 
313 {
314  assert(variable_id < getNumberOfVariables());
315  return _variable_component_offsets[variable_id + 1] -
316  _variable_component_offsets[variable_id];
317 }
318 
320 {
321  return _mesh_subsets.size();
322 }
323 
324 std::size_t LocalToGlobalIndexMap::size() const
325 {
326  return _rows.rows();
327 }
328 
330  std::size_t const mesh_item_id, const int component_id) const
331 {
332  return RowColumnIndices(_rows(mesh_item_id, component_id),
333  _columns(mesh_item_id, component_id));
334 }
335 
336 std::size_t
337 LocalToGlobalIndexMap::getNumberOfElementDOF(std::size_t const mesh_item_id) const
338 {
339  std::size_t ndof = 0;
340 
341  for (Table::Index c = 0; c < _rows.cols(); ++c)
342  {
343  ndof += _rows(mesh_item_id, c).size();
344  }
345 
346  return ndof;
347 }
348 
349 std::size_t
350 LocalToGlobalIndexMap::getNumberOfElementComponents(std::size_t const mesh_item_id) const
351 {
352  std::size_t n = 0;
353  for (Table::Index c = 0; c < _rows.cols(); ++c)
354  {
355  if (!_rows(mesh_item_id, c).empty())
356  n++;
357  }
358  return n;
359 }
360 
362  std::size_t const mesh_item_id) const
363 {
364  std::vector<int> vec;
365  for (int i = 0; i < getNumberOfVariables(); i++)
366  {
367  for (int j=0; j<getNumberOfVariableComponents(i); j++)
368  {
369  auto comp_id = getGlobalComponent(i, j);
370  if (!_rows(mesh_item_id, comp_id).empty())
371  vec.push_back(i);
372  }
373  }
374  std::sort(vec.begin(), vec.end());
375  vec.erase(std::unique(vec.begin(), vec.end()), vec.end());
376 
377  return vec;
378 }
379 
381  MeshLib::Location const& l,
382  int const variable_id,
383  int const component_id) const
384 {
385  auto const c = getGlobalComponent(variable_id, component_id);
387 }
388 
390  MeshLib::Location const& l, int const global_component_id) const
391 {
392  return _mesh_component_map.getGlobalIndex(l, global_component_id);
393 }
394 
396 std::vector<GlobalIndexType> LocalToGlobalIndexMap::getGlobalIndices(
397  const MeshLib::Location& l) const
398 {
400 }
401 
403 std::vector<GlobalIndexType> const& LocalToGlobalIndexMap::getGhostIndices()
404  const
405 {
407 }
408 
412  MeshLib::Location const& l, std::size_t const comp_id,
413  std::size_t const range_begin, std::size_t const range_end) const
414 {
415  return _mesh_component_map.getLocalIndex(l, comp_id, range_begin,
416  range_end);
417 }
418 
420  int const variable_id, int const component_id) const
421 {
422  return getMeshSubset(getGlobalComponent(variable_id, component_id));
423 }
424 
426  int const global_component_id) const
427 {
428  return _mesh_subsets[global_component_id];
429 }
430 
431 #ifndef NDEBUG
432 std::ostream& operator<<(std::ostream& os, LocalToGlobalIndexMap const& map)
433 {
434  std::size_t const max_lines = 10;
435  std::size_t lines_printed = 0;
436 
437  os << "Rows of the local to global index map; " << map._rows.size()
438  << " rows\n";
439  for (std::size_t e=0; e<map.size(); ++e)
440  {
441  os << "== e " << e << " ==\n";
442  for (int c = 0; c < map.getNumberOfComponents(); ++c)
443  {
444  auto const& line = map._rows(e, c);
445 
446  os << "c" << c << " { ";
447  std::copy(line.cbegin(), line.cend(),
448  std::ostream_iterator<std::size_t>(os, " "));
449  os << " }\n";
450  }
451 
452  if (lines_printed++ > max_lines)
453  {
454  os << "...\n";
455  break;
456  }
457  }
458  lines_printed = 0;
459 
460  os << "Mesh component map:\n" << map._mesh_component_map;
461  return os;
462 }
463 #endif // NDEBUG
464 
465 } // 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:71
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