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