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 
18 #ifdef USE_PETSC
20 #endif
21 
22 namespace NumLib
23 {
24 using namespace detail;
25 
27  std::numeric_limits<GlobalIndexType>::max();
28 
29 #ifdef USE_PETSC
31  std::vector<MeshLib::MeshSubset> const& components, ComponentOrder order)
32 {
33  // Use PETSc with single thread
34  const MeshLib::NodePartitionedMesh& partitioned_mesh =
35  static_cast<const MeshLib::NodePartitionedMesh&>(
36  components[0].getMesh());
37  if (partitioned_mesh.isForSingleThread())
38  {
39  createSerialMeshComponentMap(components, order);
40  return;
41  }
42  else
43  {
44  createParallelMeshComponentMap(components, order);
45  }
46 }
47 #else
49  std::vector<MeshLib::MeshSubset> const& components, ComponentOrder order)
50 {
51  createSerialMeshComponentMap(components, order);
52 }
53 #endif // end of USE_PETSC
54 
56  std::vector<MeshLib::MeshSubset> const& bulk_mesh_subsets,
57  MeshLib::MeshSubset const& new_mesh_subset,
58  std::vector<int> const& new_global_component_ids) const
59 {
60  { // Testing first an assumption met later in the code that the meshes for
61  // all the bulk_mesh_subsets are equal.
62  auto const first_mismatch =
63  std::adjacent_find(begin(bulk_mesh_subsets), end(bulk_mesh_subsets),
64  [](auto const& a, auto const& b)
65  { return a.getMeshID() != b.getMeshID(); });
66  if (first_mismatch != end(bulk_mesh_subsets))
67  {
68  OGS_FATAL(
69  "Assumption in the MeshComponentMap violated. Expecting all of "
70  "mesh ids to be the same, but it is not true for "
71  "the mesh '{:s}' with id {:d}.",
72  first_mismatch->getMesh().getName(),
73  first_mismatch->getMeshID());
74  }
75  }
76 
77  // Mapping of the nodes in the new_mesh_subset to the bulk mesh nodes
78  auto const& new_mesh_properties = new_mesh_subset.getMesh().getProperties();
79  if (!new_mesh_properties.template existsPropertyVector<std::size_t>(
80  "bulk_node_ids"))
81  {
82  OGS_FATAL(
83  "Bulk node ids map expected in the construction of the mesh "
84  "subset.");
85  }
86  auto const& bulk_node_ids_map =
87  *new_mesh_properties.template getPropertyVector<std::size_t>(
88  "bulk_node_ids", MeshLib::MeshItemType::Node, 1);
89 
90  // New dictionary for the subset.
91  ComponentGlobalIndexDict subset_dict;
92 
93  std::size_t const new_mesh_id = new_mesh_subset.getMeshID();
94  // Lookup the locations in the current mesh component map and
95  // insert the full lines into the new subset dictionary.
96  for (auto* const node : new_mesh_subset.getNodes())
97  {
98  auto const node_id = node->getID();
99  bool const is_base_node = MeshLib::isBaseNode(
100  *node,
101  new_mesh_subset.getMesh().getElementsConnectedToNode(node_id));
102 
103  MeshLib::Location const new_location{
104  new_mesh_id, MeshLib::MeshItemType::Node, node_id};
105 
106  // Assuming the meshes for all the bulk_mesh_subsets are equal.
107  MeshLib::Location const bulk_location{
108  bulk_mesh_subsets.front().getMeshID(), MeshLib::MeshItemType::Node,
109  bulk_node_ids_map[node_id]};
110 
111  for (auto component_id : new_global_component_ids)
112  {
113  auto const global_index =
114  getGlobalIndex(bulk_location, component_id);
115  if (global_index == nop)
116  {
117  if (is_base_node)
118  {
119  OGS_FATAL(
120  "Could not find a global index for global component "
121  "{:d} for the mesh '{:s}', node {:d}, in the "
122  "corresponding bulk mesh '{:s}' and node {:d}. This "
123  "happens because the boundary mesh is larger then the "
124  "definition region of the bulk component, usually "
125  "because the geometry for the boundary condition is "
126  "too large.",
127  component_id,
128  new_mesh_subset.getMesh().getName(),
129  node_id,
130  bulk_mesh_subsets.front().getMesh().getName(),
131  bulk_node_ids_map[node_id]);
132  }
133  continue;
134  }
135  subset_dict.insert({new_location, component_id, global_index});
136  }
137  }
138 
139  return MeshComponentMap(subset_dict);
140 }
141 
143 {
144  GlobalIndexType global_index = offset;
145 
146  auto& m = _dict.get<ByLocation>(); // view as sorted by mesh item
147  for (auto itr_mesh_item = m.begin(); itr_mesh_item != m.end();
148  ++itr_mesh_item)
149  {
150  Line pos = *itr_mesh_item;
151  pos.global_index = global_index++;
152  m.replace(itr_mesh_item, pos);
153  }
154 }
155 
156 std::vector<int> MeshComponentMap::getComponentIDs(const Location& l) const
157 {
158  auto const& m = _dict.get<ByLocation>();
159  auto const p = m.equal_range(Line(l));
160  std::vector<int> vec_compID;
161  for (auto itr = p.first; itr != p.second; ++itr)
162  {
163  vec_compID.push_back(itr->comp_id);
164  }
165  return vec_compID;
166 }
167 
168 Line MeshComponentMap::getLine(Location const& l, int const comp_id) const
169 {
170  auto const& m = _dict.get<ByLocationAndComponent>();
171  auto const itr = m.find(Line(l, comp_id));
172  assert(itr != m.end()); // The line must exist in the current dictionary.
173  return *itr;
174 }
175 
177  int const comp_id) const
178 {
179  auto const& m = _dict.get<ByLocationAndComponent>();
180  auto const itr = m.find(Line(l, comp_id));
181  return itr != m.end() ? itr->global_index : nop;
182 }
183 
184 std::vector<GlobalIndexType> MeshComponentMap::getGlobalIndices(
185  const Location& l) const
186 {
187  auto const& m = _dict.get<ByLocation>();
188  auto const p = m.equal_range(Line(l));
189  std::vector<GlobalIndexType> global_indices;
190  for (auto itr = p.first; itr != p.second; ++itr)
191  {
192  global_indices.push_back(itr->global_index);
193  }
194  return global_indices;
195 }
196 
197 std::vector<GlobalIndexType> MeshComponentMap::getGlobalIndicesByLocation(
198  std::vector<Location> const& ls) const
199 {
200  // Create vector of global indices sorted by location containing all
201  // locations given in ls parameter.
202 
203  std::vector<GlobalIndexType> global_indices;
204  global_indices.reserve(ls.size());
205 
206  auto const& m = _dict.get<ByLocation>();
207  for (const auto& l : ls)
208  {
209  auto const p = m.equal_range(Line(l));
210  for (auto itr = p.first; itr != p.second; ++itr)
211  {
212  global_indices.push_back(itr->global_index);
213  }
214  }
215 
216  return global_indices;
217 }
218 
220  std::vector<Location> const& ls) const
221 {
222  // vector of (Component, global Index) pairs.
223  using CIPair = std::pair<int, GlobalIndexType>;
224  std::vector<CIPair> pairs;
225  pairs.reserve(ls.size());
226 
227  // Create a sub dictionary containing all lines with location from ls.
228  auto const& m = _dict.get<ByLocation>();
229  for (const auto& l : ls)
230  {
231  auto const p = m.equal_range(Line(l));
232  for (auto itr = p.first; itr != p.second; ++itr)
233  {
234  pairs.emplace_back(itr->comp_id, itr->global_index);
235  }
236  }
237 
238  auto CIPairLess = [](CIPair const& a, CIPair const& b)
239  { return a.first < b.first; };
240 
241  // Create vector of global indices from sub dictionary sorting by component.
242  if (!std::is_sorted(pairs.begin(), pairs.end(), CIPairLess))
243  {
244  std::stable_sort(pairs.begin(), pairs.end(), CIPairLess);
245  }
246 
247  std::vector<GlobalIndexType> global_indices;
248  global_indices.reserve(pairs.size());
249  transform(cbegin(pairs), cend(pairs), back_inserter(global_indices),
250  [&](const auto& pair) { return pair.second; });
251 
252  return global_indices;
253 }
254 
256  Location const& l,
257  int const comp_id,
258  std::size_t const range_begin,
259  std::size_t const range_end) const
260 {
261  GlobalIndexType const global_index = getGlobalIndex(l, comp_id);
262 #ifndef USE_PETSC
263  (void)range_begin;
264  (void)range_end;
265  return global_index;
266 #else
267  if (global_index >= 0) // non-ghost location.
268  return global_index - range_begin;
269 
270  //
271  // For a ghost location look up the global index in ghost indices.
272  //
273 
274  // A special case for a ghost location with global index equal to the size
275  // of the local vector:
276  GlobalIndexType const real_global_index =
277  (-global_index == static_cast<GlobalIndexType>(_num_global_dof))
278  ? 0
279  : -global_index;
280 
281  // TODO Find in ghost indices is O(n^2/2) for n being the length of
282  // _ghosts_indices. Providing an inverted table would be faster.
283  auto const ghost_index_it = std::find(
284  _ghosts_indices.begin(), _ghosts_indices.end(), real_global_index);
285  if (ghost_index_it == _ghosts_indices.end())
286  {
287  OGS_FATAL("index {:d} not found in ghost_indices", real_global_index);
288  }
289 
290  // Using std::distance on a std::vector is O(1). As long as _ghost_indices
291  // remains of std::vector type, this shall be fast.
292  return range_end - range_begin +
293  std::distance(_ghosts_indices.begin(), ghost_index_it);
294 
295 #endif
296 }
297 
299  std::vector<MeshLib::MeshSubset> const& components, ComponentOrder order)
300 {
301  // construct dict (and here we number global_index by component type)
302  GlobalIndexType global_index = 0;
303  int comp_id = 0;
304  for (auto const& c : components)
305  {
306  std::size_t const mesh_id = c.getMeshID();
307  // mesh items are ordered first by node, cell, ....
308  for (std::size_t j = 0; j < c.getNumberOfNodes(); j++)
309  {
310  _dict.insert(Line(
311  Location(mesh_id, MeshLib::MeshItemType::Node, c.getNodeID(j)),
312  comp_id, global_index++));
313  }
314  comp_id++;
315  }
316  _num_local_dof = _dict.size();
317 
318  if (order == ComponentOrder::BY_LOCATION)
319  {
320  renumberByLocation();
321  }
322 }
323 
324 #ifdef USE_PETSC
326  std::vector<MeshLib::MeshSubset> const& components, ComponentOrder order)
327 {
328  if (order != ComponentOrder::BY_LOCATION)
329  {
330  // Not allowed in parallel case since this is not suitable to
331  // arrange non ghost entries of a partition within
332  // a rank in the parallel computing.
333  OGS_FATAL(
334  "Global index in the system of equations can only be numbered by "
335  "the order type of ComponentOrder::BY_LOCATION");
336  }
337 
338  // get number of unknowns
339  GlobalIndexType num_unknowns = 0;
340  for (auto const& c : components)
341  {
342  const MeshLib::NodePartitionedMesh& partitioned_mesh =
343  static_cast<const MeshLib::NodePartitionedMesh&>(c.getMesh());
344  num_unknowns += partitioned_mesh.getNumberOfGlobalNodes();
345  }
346 
347  // construct dict (and here we number global_index by component type)
348  int comp_id = 0;
349  _num_global_dof = 0;
350  _num_local_dof = 0;
351  for (auto const& c : components)
352  {
353  assert(dynamic_cast<MeshLib::NodePartitionedMesh const*>(
354  &c.getMesh()) != nullptr);
355  std::size_t const mesh_id = c.getMeshID();
356  const MeshLib::NodePartitionedMesh& partitioned_mesh =
357  static_cast<const MeshLib::NodePartitionedMesh&>(c.getMesh());
358 
359  // mesh items are ordered first by node, cell, ....
360  for (std::size_t j = 0; j < c.getNumberOfNodes(); j++)
361  {
362  GlobalIndexType global_index = static_cast<GlobalIndexType>(
363  components.size() * partitioned_mesh.getGlobalNodeID(j) +
364  comp_id);
365  const bool is_ghost = partitioned_mesh.isGhostNode(
366  partitioned_mesh.getNode(j)->getID());
367  if (is_ghost)
368  {
369  _ghosts_indices.push_back(global_index);
370  global_index = -global_index;
371  // If the ghost entry has an index of 0,
372  // its index is set to the negative value of unknowns.
373  if (global_index == 0)
374  {
375  global_index = -num_unknowns;
376  }
377  }
378  else
379  {
380  _num_local_dof++;
381  }
382 
383  _dict.insert(Line(Location(mesh_id, MeshLib::MeshItemType::Node, j),
384  comp_id, global_index));
385  }
386 
387  _num_global_dof += partitioned_mesh.getNumberOfGlobalNodes();
388  comp_id++;
389  }
390 }
391 #endif
392 
393 } // 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...
std::size_t getID() const
Definition: Point3dWithID.h:62
A subset of nodes on a single mesh.
Definition: MeshSubset.h:27
Mesh const & getMesh() const
Definition: MeshSubset.h:106
std::size_t getMeshID() const
return this mesh ID
Definition: MeshSubset.h:71
std::vector< Node * > const & getNodes() const
Definition: MeshSubset.h:101
const std::string getName() const
Get name of the mesh.
Definition: Mesh.h:92
Properties & getProperties()
Definition: Mesh.h:123
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
Definition: Mesh.cpp:232
const Node * getNode(std::size_t idx) const
Get the node with the given index.
Definition: Mesh.h:74
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 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.
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:31
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
ComponentOrder
Ordering of components in global matrix/vector.
@ BY_LOCATION
Ordering data by spatial location.