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