OGS
MeshComponentMap.cpp
Go to the documentation of this file.
1 
13 #include "MeshComponentMap.h"
14 
15 #include "BaseLib/Error.h"
16 #include "MeshLib/MeshSubset.h"
17 #include "MeshLib/Node.h"
18 
19 #ifdef USE_PETSC
21 #endif
22 
23 namespace NumLib
24 {
25 using namespace detail;
26 
28  std::numeric_limits<GlobalIndexType>::max();
29 
30 #ifdef USE_PETSC
32  std::vector<MeshLib::MeshSubset> const& components, ComponentOrder order)
33 {
34  // Use PETSc with single thread
35  const MeshLib::NodePartitionedMesh& partitioned_mesh =
36  static_cast<const MeshLib::NodePartitionedMesh&>(
37  components[0].getMesh());
38  if (partitioned_mesh.isForSingleThread())
39  {
40  createSerialMeshComponentMap(components, order);
41  return;
42  }
43  else
44  {
45  createParallelMeshComponentMap(components, order);
46  }
47 }
48 #else
50  std::vector<MeshLib::MeshSubset> const& components, ComponentOrder order)
51 {
52  createSerialMeshComponentMap(components, order);
53 }
54 #endif // end of USE_PETSC
55 
57  std::vector<MeshLib::MeshSubset> const& bulk_mesh_subsets,
58  MeshLib::MeshSubset const& new_mesh_subset,
59  std::vector<int> const& new_global_component_ids) const
60 {
61  { // Testing first an assumption met later in the code that the meshes for
62  // all the bulk_mesh_subsets are equal.
63  auto const first_mismatch =
64  std::adjacent_find(begin(bulk_mesh_subsets), end(bulk_mesh_subsets),
65  [](auto const& a, auto const& b)
66  { return a.getMeshID() != b.getMeshID(); });
67  if (first_mismatch != end(bulk_mesh_subsets))
68  {
69  OGS_FATAL(
70  "Assumption in the MeshComponentMap violated. Expecting all of "
71  "mesh ids to be the same, but it is not true for "
72  "the mesh '{:s}' with id {:d}.",
73  first_mismatch->getMesh().getName(),
74  first_mismatch->getMeshID());
75  }
76  }
77 
78  // Mapping of the nodes in the new_mesh_subset to the bulk mesh nodes
79  auto const& new_mesh_properties = new_mesh_subset.getMesh().getProperties();
80  if (!new_mesh_properties.template existsPropertyVector<std::size_t>(
81  "bulk_node_ids"))
82  {
83  OGS_FATAL(
84  "Bulk node ids map expected in the construction of the mesh "
85  "subset.");
86  }
87  auto const& bulk_node_ids_map =
88  *new_mesh_properties.template getPropertyVector<std::size_t>(
89  "bulk_node_ids", MeshLib::MeshItemType::Node, 1);
90 
91  // New dictionary for the subset.
92  ComponentGlobalIndexDict subset_dict;
93 
94  std::size_t const new_mesh_id = new_mesh_subset.getMeshID();
95  // Lookup the locations in the current mesh component map and
96  // insert the full lines into the new subset dictionary.
97  for (auto* const node : new_mesh_subset.getNodes())
98  {
99  auto const node_id = node->getID();
100  bool const is_base_node = MeshLib::isBaseNode(
101  *node,
102  new_mesh_subset.getMesh().getElementsConnectedToNode(node_id));
103 
104  MeshLib::Location const new_location{
105  new_mesh_id, MeshLib::MeshItemType::Node, node_id};
106 
107  // Assuming the meshes for all the bulk_mesh_subsets are equal.
108  MeshLib::Location const bulk_location{
109  bulk_mesh_subsets.front().getMeshID(), MeshLib::MeshItemType::Node,
110  bulk_node_ids_map[node_id]};
111 
112  for (auto component_id : new_global_component_ids)
113  {
114  auto const global_index =
115  getGlobalIndex(bulk_location, component_id);
116  if (global_index == nop)
117  {
118  if (is_base_node)
119  {
120  OGS_FATAL(
121  "Could not find a global index for global component "
122  "{:d} for the mesh '{:s}', node {:d}, in the "
123  "corresponding bulk mesh '{:s}' and node {:d}. This "
124  "happens because the boundary mesh is larger then the "
125  "definition region of the bulk component, usually "
126  "because the geometry for the boundary condition is "
127  "too large.",
128  component_id,
129  new_mesh_subset.getMesh().getName(),
130  node_id,
131  bulk_mesh_subsets.front().getMesh().getName(),
132  bulk_node_ids_map[node_id]);
133  }
134  continue;
135  }
136  subset_dict.insert({new_location, component_id, global_index});
137  }
138  }
139 
140  return MeshComponentMap(subset_dict);
141 }
142 
144 {
145  GlobalIndexType global_index = offset;
146 
147  auto& m = _dict.get<ByLocation>(); // view as sorted by mesh item
148  for (auto itr_mesh_item = m.begin(); itr_mesh_item != m.end();
149  ++itr_mesh_item)
150  {
151  Line pos = *itr_mesh_item;
152  pos.global_index = global_index++;
153  m.replace(itr_mesh_item, pos);
154  }
155 }
156 
157 std::vector<int> MeshComponentMap::getComponentIDs(const Location& l) const
158 {
159  auto const& m = _dict.get<ByLocation>();
160  auto const p = m.equal_range(Line(l));
161  std::vector<int> vec_compID;
162  for (auto itr = p.first; itr != p.second; ++itr)
163  {
164  vec_compID.push_back(itr->comp_id);
165  }
166  return vec_compID;
167 }
168 
169 Line MeshComponentMap::getLine(Location const& l, int const comp_id) const
170 {
171  auto const& m = _dict.get<ByLocationAndComponent>();
172  auto const itr = m.find(Line(l, comp_id));
173  assert(itr != m.end()); // The line must exist in the current dictionary.
174  return *itr;
175 }
176 
178  int const comp_id) const
179 {
180  auto const& m = _dict.get<ByLocationAndComponent>();
181  auto const itr = m.find(Line(l, comp_id));
182  return itr != m.end() ? itr->global_index : nop;
183 }
184 
185 std::vector<GlobalIndexType> MeshComponentMap::getGlobalIndices(
186  const Location& l) const
187 {
188  auto const& m = _dict.get<ByLocation>();
189  auto const p = m.equal_range(Line(l));
190  std::vector<GlobalIndexType> global_indices;
191  for (auto itr = p.first; itr != p.second; ++itr)
192  {
193  global_indices.push_back(itr->global_index);
194  }
195  return global_indices;
196 }
197 
198 std::vector<GlobalIndexType> MeshComponentMap::getGlobalIndicesByLocation(
199  std::vector<Location> const& ls) const
200 {
201  // Create vector of global indices sorted by location containing all
202  // locations given in ls parameter.
203 
204  std::vector<GlobalIndexType> global_indices;
205  global_indices.reserve(ls.size());
206 
207  auto const& m = _dict.get<ByLocation>();
208  for (const auto& l : ls)
209  {
210  auto const p = m.equal_range(Line(l));
211  for (auto itr = p.first; itr != p.second; ++itr)
212  {
213  global_indices.push_back(itr->global_index);
214  }
215  }
216 
217  return global_indices;
218 }
219 
221  std::vector<Location> const& ls) const
222 {
223  // vector of (Component, global Index) pairs.
224  using CIPair = std::pair<int, GlobalIndexType>;
225  std::vector<CIPair> pairs;
226  pairs.reserve(ls.size());
227 
228  // Create a sub dictionary containing all lines with location from ls.
229  auto const& m = _dict.get<ByLocation>();
230  for (const auto& l : ls)
231  {
232  auto const p = m.equal_range(Line(l));
233  for (auto itr = p.first; itr != p.second; ++itr)
234  {
235  pairs.emplace_back(itr->comp_id, itr->global_index);
236  }
237  }
238 
239  auto CIPairLess = [](CIPair const& a, CIPair const& b)
240  { return a.first < b.first; };
241 
242  // Create vector of global indices from sub dictionary sorting by component.
243  if (!std::is_sorted(pairs.begin(), pairs.end(), CIPairLess))
244  {
245  std::stable_sort(pairs.begin(), pairs.end(), CIPairLess);
246  }
247 
248  std::vector<GlobalIndexType> global_indices;
249  global_indices.reserve(pairs.size());
250  transform(cbegin(pairs), cend(pairs), back_inserter(global_indices),
251  [&](const auto& pair) { return pair.second; });
252 
253  return global_indices;
254 }
255 
257  Location const& l,
258  int const comp_id,
259  std::size_t const range_begin,
260  std::size_t const range_end) const
261 {
262  GlobalIndexType const global_index = getGlobalIndex(l, comp_id);
263 #ifndef USE_PETSC
264  (void)range_begin;
265  (void)range_end;
266  return global_index;
267 #else
268  if (global_index >= 0) // non-ghost location.
269  return global_index - range_begin;
270 
271  //
272  // For a ghost location look up the global index in ghost indices.
273  //
274 
275  // A special case for a ghost location with global index equal to the size
276  // of the local vector:
277  GlobalIndexType const real_global_index =
278  (-global_index == static_cast<GlobalIndexType>(_num_global_dof))
279  ? 0
280  : -global_index;
281 
282  // TODO Find in ghost indices is O(n^2/2) for n being the length of
283  // _ghosts_indices. Providing an inverted table would be faster.
284  auto const ghost_index_it = std::find(
285  _ghosts_indices.begin(), _ghosts_indices.end(), real_global_index);
286  if (ghost_index_it == _ghosts_indices.end())
287  {
288  OGS_FATAL("index {:d} not found in ghost_indices", real_global_index);
289  }
290 
291  // Using std::distance on a std::vector is O(1). As long as _ghost_indices
292  // remains of std::vector type, this shall be fast.
293  return range_end - range_begin +
294  std::distance(_ghosts_indices.begin(), ghost_index_it);
295 
296 #endif
297 }
298 
300  std::vector<MeshLib::MeshSubset> const& components, ComponentOrder order)
301 {
302  // construct dict (and here we number global_index by component type)
303  GlobalIndexType global_index = 0;
304  int comp_id = 0;
305  for (auto const& c : components)
306  {
307  std::size_t const mesh_id = c.getMeshID();
308  // mesh items are ordered first by node, cell, ....
309  for (std::size_t j = 0; j < c.getNumberOfNodes(); j++)
310  {
311  _dict.insert(Line(
312  Location(mesh_id, MeshLib::MeshItemType::Node, c.getNodeID(j)),
313  comp_id, global_index++));
314  }
315  comp_id++;
316  }
317  _num_local_dof = _dict.size();
318 
319  if (order == ComponentOrder::BY_LOCATION)
320  {
321  renumberByLocation();
322  }
323 }
324 
325 #ifdef USE_PETSC
326 
328  MeshLib::NodePartitionedMesh const& partitioned_mesh,
329  std::size_t const global_node_id,
330  int const number_of_components_at_base_node,
331  int const number_of_components_at_high_order_node,
332  int const component_id_at_base_node,
333  int const component_id_at_high_order_node,
334  bool const is_base_node)
335 {
336  int const partition_id = partitioned_mesh.getPartitionID(global_node_id);
337 
338  auto const n_total_active_base_nodes_before_this_rank =
339  partitioned_mesh.getNumberOfActiveBaseNodesAtRank(partition_id);
340 
341  auto const n_total_active_high_order_nodes_before_this_rank =
342  partitioned_mesh.getNumberOfActiveHighOrderNodesAtRank(partition_id);
343 
344  auto const node_id_offset =
345  n_total_active_base_nodes_before_this_rank +
346  n_total_active_high_order_nodes_before_this_rank;
347 
348  auto const index_offset = n_total_active_base_nodes_before_this_rank *
349  number_of_components_at_base_node +
350  n_total_active_high_order_nodes_before_this_rank *
351  number_of_components_at_high_order_node;
352 
353  if (is_base_node)
354  {
355  return static_cast<GlobalIndexType>(
356  index_offset + (global_node_id - node_id_offset) *
357  number_of_components_at_base_node) +
358  component_id_at_base_node;
359  }
360 
361  int const n_active_base_nodes_of_this_partition =
362  partitioned_mesh.getNumberOfActiveBaseNodesAtRank(partition_id + 1) -
363  n_total_active_base_nodes_before_this_rank;
364 
365  /*
366  The global indices of components are numbered as what depicted below by
367  assuming that the base node has three components and the high order node
368  has two components:
369 
370  Partition | 0 | 1 | ...
371  --------------------------------
372  Active nodes | Base | higher | Base | higher| ...
373  --------------------------------
374  c0 x x x x ...
375  c1 x x x x ...
376  c2 x x ...
377  */
378  return static_cast<GlobalIndexType>(
379  index_offset +
380  n_active_base_nodes_of_this_partition *
381  number_of_components_at_base_node +
382  (global_node_id - node_id_offset -
383  n_active_base_nodes_of_this_partition) *
384  number_of_components_at_high_order_node) +
385  component_id_at_high_order_node;
386 }
387 
389  std::vector<MeshLib::MeshSubset> const& components, ComponentOrder order)
390 {
391  if (order != ComponentOrder::BY_LOCATION)
392  {
393  // Not allowed in parallel case since this is not suitable to
394  // arrange non ghost entries of a partition within
395  // a rank in the parallel computing.
396  OGS_FATAL(
397  "Global index in the system of equations can only be numbered by "
398  "the order type of ComponentOrder::BY_LOCATION");
399  }
400 
401  const MeshLib::NodePartitionedMesh& partitioned_mesh =
402  static_cast<const MeshLib::NodePartitionedMesh&>(
403  components[0].getMesh());
404 
405  //
406  // get the number of unknowns and the number of components at extra node
407  //
408  GlobalIndexType num_unknowns = 0;
409  int components_at_high_order_nodes = 0;
410  for (auto const& c : components)
411  {
412  if (partitioned_mesh.getNumberOfNodes() == c.getNumberOfNodes())
413  {
414  components_at_high_order_nodes++;
415  num_unknowns += partitioned_mesh.getNumberOfGlobalNodes();
416  }
417  else
418  {
419  num_unknowns += partitioned_mesh.getNumberOfGlobalBaseNodes();
420  }
421  }
422 
423  // construct dict (and here we number global_index by component type)
424  //
425  int const n_components = components.size();
426 
427  int comp_id = 0;
428  int comp_id_at_high_order_node = 0;
429  _num_global_dof = 0;
430  _num_local_dof = 0;
431  for (auto const& c : components)
432  {
433  assert(dynamic_cast<MeshLib::NodePartitionedMesh const*>(
434  &c.getMesh()) != nullptr);
435  std::size_t const mesh_id = c.getMeshID();
436  const MeshLib::NodePartitionedMesh& partitioned_mesh =
437  static_cast<const MeshLib::NodePartitionedMesh&>(c.getMesh());
438  const auto& sub_mesh_nodes = c.getNodes();
439 
440  // mesh items are ordered first by node, cell, ....
441  for (std::size_t j = 0; j < c.getNumberOfNodes(); j++)
442  {
443  const auto node_id = c.getNodeID(j);
444  const bool is_base_node =
445  MeshLib::isBaseNode(*sub_mesh_nodes[j],
446  partitioned_mesh.getElementsConnectedToNode(
447  *sub_mesh_nodes[j]));
448 
449  const auto global_node_id =
450  partitioned_mesh.getGlobalNodeID(node_id);
451  GlobalIndexType global_index =
452  (!c.useTaylorHoodElements())
453  ? static_cast<GlobalIndexType>(
454  n_components * global_node_id + comp_id)
456  partitioned_mesh, global_node_id, n_components,
457  components_at_high_order_nodes, comp_id,
458  comp_id_at_high_order_node, is_base_node);
459  const bool is_ghost = partitioned_mesh.isGhostNode(node_id);
460  if (is_ghost)
461  {
462  _ghosts_indices.push_back(global_index);
463  global_index = -global_index;
464  // If the ghost entry has an index of 0,
465  // its index is set to the negative value of unknowns.
466  if (global_index == 0)
467  {
468  global_index = -num_unknowns;
469  }
470  }
471  else
472  {
473  _num_local_dof++;
474  }
475 
476  _dict.insert(
477  Line(Location(mesh_id, MeshLib::MeshItemType::Node, node_id),
478  comp_id, global_index));
479  }
480 
481  bool const use_whole_nodes =
482  (partitioned_mesh.getNumberOfNodes() == c.getNumberOfNodes());
483  if (use_whole_nodes)
484  {
485  _num_global_dof += partitioned_mesh.getNumberOfGlobalNodes();
486  comp_id_at_high_order_node++;
487  }
488  else
489  {
490  _num_global_dof += partitioned_mesh.getNumberOfGlobalBaseNodes();
491  }
492 
493  comp_id++;
494  }
495 }
496 #endif
497 
498 } // namespace NumLib
#define OGS_FATAL(...)
Definition: Error.h:26
GlobalMatrix::IndexType GlobalIndexType
Definition of mesh class for partitioned mesh (by node) for parallel computing within the framework o...
Definition of the Node class.
A subset of nodes on a single mesh.
Definition: MeshSubset.h:26
Mesh const & getMesh() const
Definition: MeshSubset.h:104
std::size_t getMeshID() const
return this mesh ID
Definition: MeshSubset.h:76
std::vector< Node * > const & getNodes() const
Definition: MeshSubset.h:102
const std::string getName() const
Get name of the mesh.
Definition: Mesh.h:92
Properties & getProperties()
Definition: Mesh.h:123
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition: Mesh.h:89
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
Definition: Mesh.cpp:232
std::size_t getPartitionID(const std::size_t global_node_id) const
bool isGhostNode(const std::size_t node_id) const
Check whether a node with ID of node_id is a ghost node.
std::size_t getNumberOfActiveBaseNodesAtRank(int const partition_id) const
std::size_t getGlobalNodeID(const std::size_t node_id) const
Get the global node ID of a node with its local ID.
std::size_t getNumberOfGlobalNodes() const
Get the number of all nodes of the global mesh.
std::size_t getNumberOfGlobalBaseNodes() const
Get the number of nodes of the global mesh for linear elements.
std::size_t getNumberOfActiveHighOrderNodesAtRank(int const partition_id) const
Multidirectional mapping between mesh entities and degrees of freedom.
void createParallelMeshComponentMap(std::vector< MeshLib::MeshSubset > const &components, ComponentOrder order)
std::vector< int > getComponentIDs(const Location &l) const
MeshComponentMap(std::vector< MeshLib::MeshSubset > const &components, ComponentOrder order)
GlobalIndexType getGlobalIndex(Location const &l, int const comp_id) const
detail::Line getLine(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
void renumberByLocation(GlobalIndexType offset=0)
static NUMLIB_EXPORT GlobalIndexType const nop
std::vector< GlobalIndexType > getGlobalIndicesByComponent(const std::vector< Location > &ls) const
void createSerialMeshComponentMap(std::vector< MeshLib::MeshSubset > const &components, ComponentOrder order)
std::vector< GlobalIndexType > getGlobalIndicesByLocation(const std::vector< Location > &ls) const
Location
Definition: Polyline.h:33
TemplateElement< MeshLib::LineRule2 > Line
Definition: Line.h:25
bool isBaseNode(Node const &node, std::vector< Element const * > const &elements_connected_to_node)
Definition: Mesh.cpp:370
boost::multi_index::multi_index_container< Line, boost::multi_index::indexed_by< boost::multi_index::ordered_unique< boost::multi_index::tag< ByLocationAndComponent >, boost::multi_index::identity< Line >, LineByLocationAndComponentComparator >, boost::multi_index::ordered_non_unique< boost::multi_index::tag< ByLocation >, boost::multi_index::identity< Line >, LineByLocationComparator >, boost::multi_index::ordered_non_unique< boost::multi_index::tag< ByComponent >, boost::multi_index::member< Line, int, &Line::comp_id > >, boost::multi_index::ordered_non_unique< boost::multi_index::tag< ByGlobalIndex >, boost::multi_index::member< Line, GlobalIndexType, &Line::global_index > >> > ComponentGlobalIndexDict
GlobalIndexType getGlobalIndexWithTaylorHoodElement(MeshLib::NodePartitionedMesh const &partitioned_mesh, std::size_t const global_node_id, int const number_of_components_at_base_node, int const number_of_components_at_high_order_node, int const component_id_at_base_node, int const component_id_at_high_order_node, bool const is_base_node)
ComponentOrder
Ordering of components in global matrix/vector.
@ BY_LOCATION
Ordering data by spatial location.