OGS 6.1.0-1721-g6382411ad
IdentifySubdomainMesh.cpp
Go to the documentation of this file.
1 
10 #include <map>
11 #include <vector>
12 
14 #include "MeshLib/Mesh.h"
15 #include "MeshLib/Node.h"
16 
17 #include "MeshNodeSearcher.h"
18 
19 namespace
20 {
23 std::vector<std::size_t> identifySubdomainMeshNodes(
24  MeshLib::Mesh const& subdomain_mesh,
25  MeshGeoToolsLib::MeshNodeSearcher const& mesh_node_searcher)
26 {
27  // Convert nodes pointers needed for the mesh_node_searcher algorithm.
28  auto const& nodes = subdomain_mesh.getNodes();
29  std::vector<MathLib::Point3dWithID*> subdomain_points{begin(nodes),
30  end(nodes)};
31 
32  auto const& bulk_node_ids =
33  mesh_node_searcher.getMeshNodeIDs(subdomain_points);
34 
35  if (bulk_node_ids.size() != subdomain_mesh.getNumberOfNodes())
36  {
37  OGS_FATAL(
38  "Expected to find exactly one node in the bulk mesh for each node "
39  "of the subdomain; Found %d nodes in the bulk mesh out of %d nodes "
40  "in the subdomain.",
41  bulk_node_ids.size(), subdomain_mesh.getNumberOfNodes());
42  }
43 
44  return bulk_node_ids;
45 }
46 
50 std::vector<std::size_t> findElementsInMesh(
51  MeshLib::Mesh const& mesh, std::vector<std::size_t> const& node_ids)
52 {
53  //
54  // Collect all element ids for all nodes.
55  //
56  std::vector<std::size_t> common_element_ids;
57  // Every node is connected to at least one element.
58  auto const nnodes = node_ids.size();
59  common_element_ids.reserve(nnodes);
60 
61  for (auto const node_id : node_ids)
62  {
63  auto const& connected_elements = mesh.getNode(node_id)->getElements();
64  std::transform(
65  begin(connected_elements), end(connected_elements),
66  back_inserter(common_element_ids),
67  [](MeshLib::Element const* const e) { return e->getID(); });
68  }
69 
70  //
71  // Count how often an element is shared by all nodes.
72  //
73  std::map<std::size_t, int> element_counts;
74  for (auto const element_id : common_element_ids)
75  {
76  element_counts[element_id]++;
77  }
78 
79  //
80  // Elements which are shared by as many nodes as the input nodes are the
81  // desired elements.
82  //
83  std::vector<std::size_t> element_ids;
84  for (auto const& pair : element_counts)
85  {
86  if (pair.second == static_cast<int>(nnodes))
87  {
88  element_ids.push_back(pair.first);
89  }
90  }
91 
92  return element_ids;
93 }
94 
98 std::vector<std::vector<std::size_t>> identifySubdomainMeshElements(
99  MeshLib::Mesh const& subdomain_mesh, MeshLib::Mesh const& bulk_mesh)
100 {
101  auto& properties = subdomain_mesh.getProperties();
102  auto const& bulk_node_ids = *properties.getPropertyVector<std::size_t>(
103  "bulk_node_ids", MeshLib::MeshItemType::Node, 1);
104 
105  // Allocate space for all elements for random insertion.
106  std::vector<std::vector<std::size_t>> bulk_element_ids_map(
107  subdomain_mesh.getNumberOfElements());
108 
109  for (auto* const e : subdomain_mesh.getElements())
110  {
111  std::vector<std::size_t> element_node_ids(e->getNumberOfBaseNodes());
112  for (unsigned n = 0; n < e->getNumberOfBaseNodes(); ++n)
113  {
114  element_node_ids[n] = e->getNodeIndex(n);
115  }
116  std::vector<std::size_t> element_node_ids_bulk(
117  e->getNumberOfBaseNodes());
118  std::transform(begin(element_node_ids), end(element_node_ids),
119  begin(element_node_ids_bulk),
120  [&bulk_node_ids](std::size_t const id) {
121  return bulk_node_ids[id];
122  });
123  std::vector<std::size_t> bulk_element_ids =
124  findElementsInMesh(bulk_mesh, element_node_ids_bulk);
125 
126  if (bulk_element_ids.empty())
127  {
128  ERR("No element could be found for the subdomain element %d. "
129  "Corresponding bulk mesh node ids are:",
130  e->getID());
131  for (auto const i : element_node_ids_bulk)
132  ERR("\t%d", i);
133  OGS_FATAL(
134  "Expect at least one element to be found in the bulk mesh.");
135  }
136 
137  bulk_element_ids_map[e->getID()] = std::move(bulk_element_ids);
138  }
139 
140  return bulk_element_ids_map;
141 }
142 
145  MeshLib::Mesh& mesh, std::string const& property_name,
146  std::vector<std::size_t> const& values,
147  MeshLib::MeshItemType const mesh_item_type, bool const force_overwrite)
148 {
149  auto& properties = mesh.getProperties();
150  if (!properties.existsPropertyVector<std::size_t>(property_name))
151  {
152  addPropertyToMesh(mesh, property_name, mesh_item_type, 1, values);
153  return;
154  }
155 
156  //
157  // Check the existing property against new values.
158  //
159  auto& original_property =
160  *properties.getPropertyVector<std::size_t>(property_name);
161  if (std::equal(begin(original_property), end(original_property),
162  begin(values), end(values)))
163  {
164  INFO(
165  "There is already a '%s' property present in the subdomain mesh "
166  "'%s' and it is equal to the newly computed values.",
167  property_name.c_str(), mesh.getName().c_str());
168  return;
169  }
170 
171  //
172  // Property differs. Notify and update if forced.
173  //
174  WARN(
175  "There is already a '%s' property present in the subdomain mesh '%s' "
176  "and it is not equal to the newly computed values.",
177  property_name.c_str(),
178  mesh.getName().c_str());
179 
180  if (!force_overwrite)
181  {
182  OGS_FATAL("The force overwrite flag was not specified, exiting.");
183  }
184 
185  INFO("Overwriting '%s' property.", property_name.c_str());
186  original_property.resize(values.size());
187  std::copy(begin(values), end(values), begin(original_property));
188 }
189 } // namespace
190 
191 namespace MeshGeoToolsLib
192 {
194  MeshLib::Mesh const& bulk_mesh,
195  MeshNodeSearcher const& mesh_node_searcher,
196  bool const force_overwrite = false)
197 {
198  auto const& bulk_node_ids =
199  identifySubdomainMeshNodes(subdomain_mesh, mesh_node_searcher);
200 
202  subdomain_mesh, "bulk_node_ids", bulk_node_ids,
203  MeshLib::MeshItemType::Node, force_overwrite);
204 
205  auto const& bulk_element_ids =
206  identifySubdomainMeshElements(subdomain_mesh, bulk_mesh);
207 
208  // The bulk_element_ids could be of two types: one element per entry---this
209  // is the expected case for the boundary meshes; multiple elements per
210  // entry---this happens if the subdomain mesh lies inside the bulk mesh and
211  // has lower dimension.
212  // First find out the type, then add/check the CellData or FieldData.
213  if (all_of(begin(bulk_element_ids), end(bulk_element_ids),
214  [](std::vector<std::size_t> const& v) { return v.size() == 1; }))
215  {
216  // All vectors are of size 1, so the data can be flattened and
217  // stored in CellData or compared to existing CellData.
218  std::vector<std::size_t> unique_bulk_element_ids;
219  unique_bulk_element_ids.reserve(bulk_element_ids.size());
220  transform(begin(bulk_element_ids), end(bulk_element_ids),
221  back_inserter(unique_bulk_element_ids),
222  [](std::vector<std::size_t> const& v) { return v[0]; });
223 
225  subdomain_mesh, "bulk_element_ids", unique_bulk_element_ids,
226  MeshLib::MeshItemType::Cell, force_overwrite);
227  }
228  else
229  {
230  // Some of the boundary elements are connected to multiple bulk
231  // elements; Store the array in FieldData with additional CellData array
232  // for the number of elements, which also provides the offsets.
233  std::vector<std::size_t> flat_bulk_element_ids;
234  flat_bulk_element_ids.reserve(2 * bulk_element_ids.size()); // Guess.
235  std::vector<std::size_t> number_of_bulk_element_ids;
236  number_of_bulk_element_ids.reserve(bulk_element_ids.size());
237 
238  for (std::vector<std::size_t> const& v : bulk_element_ids)
239  {
240  number_of_bulk_element_ids.push_back(v.size());
241  flat_bulk_element_ids.insert(end(flat_bulk_element_ids), begin(v),
242  end(v));
243  }
244 
246  subdomain_mesh, "number_bulk_elements", number_of_bulk_element_ids,
247  MeshLib::MeshItemType::Cell, force_overwrite);
249  subdomain_mesh, "bulk_element_ids", flat_bulk_element_ids,
250  MeshLib::MeshItemType::IntegrationPoint, force_overwrite);
251  }
252 }
253 } // namespace MeshGeoToolsLib
Implementation of heuristic search length strategy.
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition: Mesh.h:99
void updateOrCheckExistingSubdomainProperty(MeshLib::Mesh &mesh, std::string const &property_name, std::vector< std::size_t > const &values, MeshLib::MeshItemType const mesh_item_type, bool const force_overwrite)
Updates or checks the existing mesh&#39;s property with the given values.
const Node * getNode(std::size_t idx) const
Get the node with the given index.
Definition: Mesh.h:84
const std::vector< Element * > & getElements() const
Get all elements the node is part of.
Definition: Node.h:65
Definition of the Node class.
std::vector< std::size_t > identifySubdomainMeshNodes(MeshLib::Mesh const &subdomain_mesh, MeshGeoToolsLib::MeshNodeSearcher const &mesh_node_searcher)
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition: Mesh.h:105
std::vector< std::vector< std::size_t > > identifySubdomainMeshElements(MeshLib::Mesh const &subdomain_mesh, MeshLib::Mesh const &bulk_mesh)
MeshLib::Properties & getProperties()
Definition: Mesh.h:134
Definition of the Mesh class.
void identifySubdomainMesh(MeshLib::Mesh &subdomain_mesh, MeshLib::Mesh const &bulk_mesh, MeshNodeSearcher const &mesh_node_searcher, bool const force_overwrite=false)
std::vector< std::size_t > findElementsInMesh(MeshLib::Mesh const &mesh, std::vector< std::size_t > const &node_ids)
const std::string getName() const
Get name of the mesh.
Definition: Mesh.h:102
MeshItemType
Definition: Location.h:21
std::vector< std::size_t > getMeshNodeIDs(GeoLib::GeoObject const &geoObj) const
PropertyVector< T > const * getPropertyVector(std::string const &name) const
Definition: Properties.h:119
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:108
void addPropertyToMesh(MeshLib::Mesh &mesh, std::string const &name, MeshLib::MeshItemType item_type, std::size_t number_of_components, std::vector< T > const &values)
Definition: Mesh.h:209
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:90
std::size_t getNumberOfElements() const
Get the number of elements.
Definition: Mesh.h:96
#define OGS_FATAL(fmt,...)
Definition: Error.h:71
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:36
Definition of the Element class.