OGS 6.1.0-1721-g6382411ad
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  _dict.insert(Line(
106  Location(mesh_id, MeshLib::MeshItemType::Node, c.getNodeID(j)),
107  comp_id, global_index++));
108  comp_id++;
109  }
110  _num_local_dof = _dict.size();
111 
112  if (order == ComponentOrder::BY_LOCATION)
113  renumberByLocation();
114 }
115 #endif // end of USE_PETSC
116 
118  std::vector<MeshLib::MeshSubset> const& bulk_mesh_subsets,
119  MeshLib::MeshSubset const& new_mesh_subset,
120  std::vector<int> const& new_global_component_ids) const
121 {
122  { // Testing first an assumption met later in the code that the meshes for
123  // the all bulk_mesh_subsets are equal.
124  auto const first_mismatch =
125  std::adjacent_find(begin(bulk_mesh_subsets), end(bulk_mesh_subsets),
126  [](auto const& a, auto const& b) {
127  return a.getMeshID() != b.getMeshID();
128  });
129  if (first_mismatch != end(bulk_mesh_subsets))
130  {
131  OGS_FATAL(
132  "Assumption in the MeshComponentMap violated. Expecting "
133  "all of mesh ids to be the same, but it is not true for "
134  "the mesh '%s' with id %d.",
135  first_mismatch->getMesh().getName().c_str(),
136  first_mismatch->getMeshID());
137  }
138  }
139 
140  // Mapping of the nodes in the new_mesh_subset to the bulk mesh nodes
141  auto const& new_mesh_properties = new_mesh_subset.getMesh().getProperties();
142  if (!new_mesh_properties.template existsPropertyVector<std::size_t>(
143  "bulk_node_ids"))
144  {
145  OGS_FATAL(
146  "Bulk node ids map expected in the construction of the mesh "
147  "subset.");
148  }
149  auto const& bulk_node_ids_map =
150  *new_mesh_properties.template getPropertyVector<std::size_t>(
151  "bulk_node_ids", MeshLib::MeshItemType::Node, 1);
152 
153  // New dictionary for the subset.
154  ComponentGlobalIndexDict subset_dict;
155 
156  std::size_t const new_mesh_id = new_mesh_subset.getMeshID();
157  // Lookup the locations in the current mesh component map and
158  // insert the full lines into the new subset dictionary.
159  for (auto* const node : new_mesh_subset.getNodes())
160  {
161  auto const node_id = node->getID();
162  bool const is_base_node = isBaseNode(*node);
163 
164  MeshLib::Location const new_location{
165  new_mesh_id, MeshLib::MeshItemType::Node, node_id};
166 
167  // Assuming the meshes for the all bulk_mesh_subsets are equal.
168  MeshLib::Location const bulk_location{
169  bulk_mesh_subsets.front().getMeshID(), MeshLib::MeshItemType::Node,
170  bulk_node_ids_map[node_id]};
171 
172  for (auto component_id : new_global_component_ids)
173  {
174  auto const global_index =
175  getGlobalIndex(bulk_location, component_id);
176  if (global_index == nop)
177  {
178  if (is_base_node)
179  {
180  OGS_FATAL(
181  "Could not find a global index for global component %d "
182  "for the mesh '%s', node %d, in the corresponding bulk "
183  "mesh '%s' and node %d. This happens because the "
184  "boundary mesh is larger then the definition region of "
185  "the bulk component, usually because the geometry for "
186  "the boundary condition is too large.",
187  component_id,
188  new_mesh_subset.getMesh().getName().c_str(),
189  node_id,
190  bulk_mesh_subsets.front().getMesh().getName().c_str(),
191  bulk_node_ids_map[node_id]);
192  }
193  continue;
194  }
195  subset_dict.insert({new_location, component_id, global_index});
196  }
197  }
198 
199  return MeshComponentMap(subset_dict);
200 }
201 
203 {
204  GlobalIndexType global_index = offset;
205 
206  auto &m = _dict.get<ByLocation>(); // view as sorted by mesh item
207  for (auto itr_mesh_item=m.begin(); itr_mesh_item!=m.end(); ++itr_mesh_item)
208  {
209  Line pos = *itr_mesh_item;
210  pos.global_index = global_index++;
211  m.replace(itr_mesh_item, pos);
212  }
213 }
214 
215 std::vector<int> MeshComponentMap::getComponentIDs(const Location& l) const
216 {
217  auto const &m = _dict.get<ByLocation>();
218  auto const p = m.equal_range(Line(l));
219  std::vector<int> vec_compID;
220  for (auto itr=p.first; itr!=p.second; ++itr)
221  vec_compID.push_back(itr->comp_id);
222  return vec_compID;
223 }
224 
225 Line MeshComponentMap::getLine(Location const& l, int const comp_id) const
226 {
227  auto const &m = _dict.get<ByLocationAndComponent>();
228  auto const itr = m.find(Line(l, comp_id));
229  assert(itr != m.end()); // The line must exist in the current dictionary.
230  return *itr;
231 }
232 
234  int const comp_id) const
235 {
236  auto const &m = _dict.get<ByLocationAndComponent>();
237  auto const itr = m.find(Line(l, comp_id));
238  return itr!=m.end() ? itr->global_index : nop;
239 }
240 
241 std::vector<GlobalIndexType> MeshComponentMap::getGlobalIndices(const Location &l) const
242 {
243  auto const &m = _dict.get<ByLocation>();
244  auto const p = m.equal_range(Line(l));
245  std::vector<GlobalIndexType> global_indices;
246  for (auto itr=p.first; itr!=p.second; ++itr)
247  global_indices.push_back(itr->global_index);
248  return global_indices;
249 }
250 
251 std::vector<GlobalIndexType> MeshComponentMap::getGlobalIndicesByLocation(
252  std::vector<Location> const& ls) const
253 {
254  // Create vector of global indices sorted by location containing all
255  // locations given in ls parameter.
256 
257  std::vector<GlobalIndexType> global_indices;
258  global_indices.reserve(ls.size());
259 
260  auto const &m = _dict.get<ByLocation>();
261  for (const auto& l : ls)
262  {
263  auto const p = m.equal_range(Line(l));
264  for (auto itr = p.first; itr != p.second; ++itr)
265  global_indices.push_back(itr->global_index);
266  }
267 
268  return global_indices;
269 }
270 
272  std::vector<Location> const& ls) const
273 {
274  // vector of (Component, global Index) pairs.
275  using CIPair = std::pair<int, GlobalIndexType>;
276  std::vector<CIPair> pairs;
277  pairs.reserve(ls.size());
278 
279  // Create a sub dictionary containing all lines with location from ls.
280  auto const &m = _dict.get<ByLocation>();
281  for (const auto& l : ls)
282  {
283  auto const p = m.equal_range(Line(l));
284  for (auto itr = p.first; itr != p.second; ++itr)
285  pairs.emplace_back(itr->comp_id, itr->global_index);
286  }
287 
288  auto CIPairLess = [](CIPair const& a, CIPair const& b)
289  {
290  return a.first < b.first;
291  };
292 
293  // Create vector of global indices from sub dictionary sorting by component.
294  if (!std::is_sorted(pairs.begin(), pairs.end(), CIPairLess))
295  std::stable_sort(pairs.begin(), pairs.end(), CIPairLess);
296 
297  std::vector<GlobalIndexType> global_indices;
298  global_indices.reserve(pairs.size());
299  for (const auto& pair : pairs)
300  global_indices.push_back(pair.second);
301 
302  return global_indices;
303 }
304 
306  Location const& l,
307  int const comp_id,
308  std::size_t const range_begin,
309  std::size_t const range_end) const
310 {
311  GlobalIndexType const global_index = getGlobalIndex(l, comp_id);
312 #ifndef USE_PETSC
313  (void)range_begin;
314  (void)range_end;
315  return global_index;
316 #else
317  if (global_index >= 0) // non-ghost location.
318  return global_index - range_begin;
319 
320  //
321  // For a ghost location look up the global index in ghost indices.
322  //
323 
324  // A special case for a ghost location with global index equal to the size
325  // of the local vector:
326  GlobalIndexType const real_global_index =
327  (-global_index == static_cast<GlobalIndexType>(_num_global_dof))
328  ? 0 : -global_index;
329 
330  // TODO Find in ghost indices is O(n^2/2) for n being the length of
331  // _ghosts_indices. Providing an inverted table would be faster.
332  auto const ghost_index_it = std::find(_ghosts_indices.begin(),
333  _ghosts_indices.end(),
334  real_global_index);
335  if (ghost_index_it == _ghosts_indices.end())
336  {
337  OGS_FATAL("index %d not found in ghost_indices",
338  real_global_index);
339  }
340 
341  // Using std::distance on a std::vector is O(1). As long as _ghost_indices
342  // remains of std::vector type, this shall be fast.
343  return range_end - range_begin +
344  std::distance(_ghosts_indices.begin(), ghost_index_it);
345 
346 #endif
347 }
348 
349 } // 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:63
bool isBaseNode(Node const &node)
Definition: Node.cpp:52
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:71
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