Loading [MathJax]/extensions/tex2jax.js
OGS
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
IdentifySubdomainMesh.cpp
Go to the documentation of this file.
1
10#include <range/v3/algorithm/equal.hpp>
11#include <range/v3/range/conversion.hpp>
12#include <unordered_map>
13#include <vector>
14
15#include "BaseLib/RunTime.h"
17#include "MeshLib/Mesh.h"
18#include "MeshLib/Node.h"
20#include "MeshNodeSearcher.h"
21
22namespace
23{
26std::vector<std::size_t> identifySubdomainMeshNodes(
27 MeshLib::Mesh const& subdomain_mesh,
28 MeshGeoToolsLib::MeshNodeSearcher const& mesh_node_searcher)
29{
30 // Convert nodes pointers needed for the mesh_node_searcher algorithm.
31 auto const& nodes = subdomain_mesh.getNodes();
32 std::vector<MathLib::Point3dWithID*> subdomain_points{begin(nodes),
33 end(nodes)};
34
35 auto const& bulk_node_ids =
36 mesh_node_searcher.getMeshNodeIDs(subdomain_points);
37
38 if (bulk_node_ids.size() != subdomain_mesh.getNumberOfNodes())
39 {
41 "Expected to find exactly one node in the bulk mesh for each node "
42 "of the subdomain; Found {:d} nodes in the bulk mesh out of {:d} "
43 "nodes in the subdomain.",
44 bulk_node_ids.size(), subdomain_mesh.getNumberOfNodes());
45 }
46
47 return bulk_node_ids;
48}
49
53std::vector<std::size_t> findElementsInMesh(
54 std::vector<std::size_t> const& node_ids,
55 std::vector<std::vector<std::size_t>> const& connected_element_ids_per_node)
56{
57 //
58 // Count how often an element is shared by all nodes.
59 //
60 std::unordered_map<std::size_t, int> element_counts(8);
61 for (auto const node_id : node_ids)
62 {
63 for (auto const element_id : connected_element_ids_per_node[node_id])
64 {
65 element_counts[element_id]++;
66 }
67 }
68
69 //
70 // Elements which are shared by as many nodes as the input nodes are the
71 // desired elements.
72 //
73 auto const nnodes = node_ids.size();
74 std::vector<std::size_t> element_ids;
75 for (auto const& pair : element_counts)
76 {
77 if (pair.second == static_cast<int>(nnodes))
78 {
79 element_ids.push_back(pair.first);
80 }
81 }
82
83 return element_ids;
84}
85
89std::vector<std::vector<std::size_t>> identifySubdomainMeshElements(
90 MeshLib::Mesh const& subdomain_mesh, MeshLib::Mesh const& bulk_mesh)
91{
92 auto const& bulk_node_ids = *MeshLib::bulkNodeIDs(subdomain_mesh);
93
94 // Allocate space for all elements for random insertion.
95 std::vector<std::vector<std::size_t>> bulk_element_ids_map(
96 subdomain_mesh.getNumberOfElements());
97
98 // For each node a vector of connected element ids of that node.
99 std::vector<std::vector<std::size_t>> connected_element_ids_per_node(
100 bulk_mesh.getNumberOfNodes());
101 for (auto const node_id : bulk_mesh.getNodes() | MeshLib::views::ids)
102 {
103 connected_element_ids_per_node[node_id] =
104 bulk_mesh.getElementsConnectedToNode(node_id) |
105 MeshLib::views::ids | ranges::to<std::vector>;
106 }
107
108 auto const& elements = subdomain_mesh.getElements();
109#pragma omp parallel for
110 for (std::ptrdiff_t j = 0; j < std::ssize(elements); ++j)
111 {
112 auto* const e = elements[j];
113 std::vector<std::size_t> element_node_ids(e->getNumberOfBaseNodes());
114 for (unsigned n = 0; n < e->getNumberOfBaseNodes(); ++n)
115 {
116 element_node_ids[n] = MeshLib::getNodeIndex(*e, n);
117 }
118 std::vector<std::size_t> element_node_ids_bulk(
119 e->getNumberOfBaseNodes());
120 std::transform(begin(element_node_ids), end(element_node_ids),
121 begin(element_node_ids_bulk),
122 [&bulk_node_ids](std::size_t const id)
123 { return bulk_node_ids[id]; });
124
125 std::vector<std::size_t> bulk_element_ids = findElementsInMesh(
126 element_node_ids_bulk, connected_element_ids_per_node);
127
128 if (bulk_element_ids.empty())
129 {
130 ERR("No element could be found for the subdomain element {:d}. "
131 "Corresponding bulk mesh node ids are:",
132 e->getID());
133 for (auto const i : element_node_ids_bulk)
134 {
135 ERR("\t{:d}", i);
136 }
137 OGS_FATAL(
138 "Expect at least one element to be found in the bulk mesh.");
139 }
140
141 bulk_element_ids_map[e->getID()] = std::move(bulk_element_ids);
142 }
143
144 return bulk_element_ids_map;
145}
146
149 MeshLib::Mesh& mesh, std::string_view property_name,
150 std::vector<std::size_t> const& values,
151 MeshLib::MeshItemType const mesh_item_type, bool const force_overwrite)
152{
153 auto& properties = mesh.getProperties();
154 if (!properties.existsPropertyVector<std::size_t>(property_name))
155 {
156 addPropertyToMesh(mesh, property_name, mesh_item_type, 1, values);
157 return;
158 }
159
160 //
161 // Check the existing property against new values.
162 //
163 auto& original_property =
164 *properties.getPropertyVector<std::size_t>(property_name);
165 if (ranges::equal(original_property, values))
166 {
167 INFO(
168 "There is already a '{:s}' property present in the subdomain mesh "
169 "'{:s}' and it is equal to the newly computed values.",
170 property_name, mesh.getName());
171 return;
172 }
173
174 //
175 // Property differs. Notify and update if forced.
176 //
177 WARN(
178 "There is already a '{:s}' property present in the subdomain mesh "
179 "'{:s}' and it is not equal to the newly computed values.",
180 property_name,
181 mesh.getName());
182
183 if (!force_overwrite)
184 {
185 OGS_FATAL("The force overwrite flag was not specified, exiting.");
186 }
187
188 INFO("Overwriting '{:s}' property.", property_name);
189 original_property.assign(values.begin(), values.end());
190}
191} // namespace
192
193namespace MeshGeoToolsLib
194{
196 MeshLib::Mesh const& bulk_mesh,
197 MeshNodeSearcher const& mesh_node_searcher,
198 bool const force_overwrite = false)
199{
200 BaseLib::RunTime time;
201 time.start();
202 auto const& bulk_node_ids =
203 identifySubdomainMeshNodes(subdomain_mesh, mesh_node_searcher);
204 INFO("identifySubdomainMesh(): identifySubdomainMeshNodes took {:g} s",
205 time.elapsed());
206
207 updateOrCheckExistingSubdomainProperty(
209 bulk_node_ids, MeshLib::MeshItemType::Node, force_overwrite);
210
211 time.start();
212 auto const& bulk_element_ids =
213 identifySubdomainMeshElements(subdomain_mesh, bulk_mesh);
214 INFO("identifySubdomainMesh(): identifySubdomainMeshElements took {:g} s",
215 time.elapsed());
216
217 // The bulk_element_ids could be of two types: one element per entry---this
218 // is the expected case for the boundary meshes; multiple elements per
219 // entry---this happens if the subdomain mesh lies inside the bulk mesh and
220 // has lower dimension.
221 // First find out the type, then add/check the CellData or FieldData.
222 if (all_of(begin(bulk_element_ids), end(bulk_element_ids),
223 [](std::vector<std::size_t> const& v) { return v.size() == 1; }))
224 {
225 // All vectors are of size 1, so the data can be flattened and
226 // stored in CellData or compared to existing CellData.
227 std::vector<std::size_t> unique_bulk_element_ids;
228 unique_bulk_element_ids.reserve(bulk_element_ids.size());
229 transform(begin(bulk_element_ids), end(bulk_element_ids),
230 back_inserter(unique_bulk_element_ids),
231 [](std::vector<std::size_t> const& v) { return v[0]; });
232
233 updateOrCheckExistingSubdomainProperty(
234 subdomain_mesh,
236 unique_bulk_element_ids, MeshLib::MeshItemType::Cell,
237 force_overwrite);
238 }
239 else
240 {
241 // Some of the boundary elements are connected to multiple bulk
242 // elements; Store the array in FieldData with additional CellData array
243 // for the number of elements, which also provides the offsets.
244 std::vector<std::size_t> flat_bulk_element_ids;
245 flat_bulk_element_ids.reserve(2 * bulk_element_ids.size()); // Guess.
246 std::vector<std::size_t> number_of_bulk_element_ids;
247 number_of_bulk_element_ids.reserve(bulk_element_ids.size());
248
249 for (std::vector<std::size_t> const& v : bulk_element_ids)
250 {
251 number_of_bulk_element_ids.push_back(v.size());
252 flat_bulk_element_ids.insert(end(flat_bulk_element_ids), begin(v),
253 end(v));
254 }
255
256 updateOrCheckExistingSubdomainProperty(
257 subdomain_mesh, "number_bulk_elements", number_of_bulk_element_ids,
258 MeshLib::MeshItemType::Cell, force_overwrite);
259 updateOrCheckExistingSubdomainProperty(
260 subdomain_mesh,
262 flat_bulk_element_ids, MeshLib::MeshItemType::IntegrationPoint,
263 force_overwrite);
264 }
265}
266} // namespace MeshGeoToolsLib
Definition of the Element class.
#define OGS_FATAL(...)
Definition Error.h:26
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40
Definition of the Mesh class.
Definition of the Node class.
Definition of the RunTime class.
Count the running time.
Definition RunTime.h:29
double elapsed() const
Get the elapsed time in seconds.
Definition RunTime.h:42
void start()
Start the timer.
Definition RunTime.h:32
std::vector< std::size_t > getMeshNodeIDs(GeoLib::GeoObject const &geoObj) const
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition Mesh.h:108
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:111
Properties & getProperties()
Definition Mesh.h:136
const std::string getName() const
Get name of the mesh.
Definition Mesh.h:105
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition Mesh.h:102
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
Definition Mesh.cpp:256
std::size_t getNumberOfElements() const
Get the number of elements.
Definition Mesh.h:99
PropertyVector< T > const * getPropertyVector(std::string_view name) const
void identifySubdomainMesh(MeshLib::Mesh &subdomain_mesh, MeshLib::Mesh const &bulk_mesh, MeshNodeSearcher const &mesh_node_searcher, bool const force_overwrite=false)
constexpr ranges::views::view_closure ids
For an element of a range view return its id.
Definition Mesh.h:227
constexpr std::string_view getBulkIDString(MeshItemType mesh_item_type)
Definition Properties.h:185
std::size_t getNodeIndex(Element const &element, unsigned const idx)
Definition Element.cpp:219
PropertyVector< std::size_t > const * bulkNodeIDs(Mesh const &mesh)
Definition Mesh.cpp:292
std::vector< std::size_t > identifySubdomainMeshNodes(MeshLib::Mesh const &subdomain_mesh, MeshGeoToolsLib::MeshNodeSearcher const &mesh_node_searcher)
void updateOrCheckExistingSubdomainProperty(MeshLib::Mesh &mesh, std::string_view 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's property with the given values.
std::vector< std::vector< std::size_t > > identifySubdomainMeshElements(MeshLib::Mesh const &subdomain_mesh, MeshLib::Mesh const &bulk_mesh)
std::vector< std::size_t > findElementsInMesh(std::vector< std::size_t > const &node_ids, std::vector< std::vector< std::size_t > > const &connected_element_ids_per_node)