10#include <range/v3/algorithm/equal.hpp>
11#include <range/v3/range/conversion.hpp>
12#include <unordered_map>
31 auto const& nodes = subdomain_mesh.
getNodes();
32 std::vector<MathLib::Point3dWithID*> subdomain_points{begin(nodes),
35 auto const& bulk_node_ids =
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.",
54 std::vector<std::size_t>
const& node_ids,
55 std::vector<std::vector<std::size_t>>
const& connected_element_ids_per_node)
60 std::unordered_map<std::size_t, int> element_counts(8);
61 for (
auto const node_id : node_ids)
63 for (
auto const element_id : connected_element_ids_per_node[node_id])
65 element_counts[element_id]++;
73 auto const nnodes = node_ids.size();
74 std::vector<std::size_t> element_ids;
75 for (
auto const& pair : element_counts)
77 if (pair.second ==
static_cast<int>(nnodes))
79 element_ids.push_back(pair.first);
95 std::vector<std::vector<std::size_t>> bulk_element_ids_map(
99 std::vector<std::vector<std::size_t>> connected_element_ids_per_node(
103 connected_element_ids_per_node[node_id] =
108 auto const& elements = subdomain_mesh.
getElements();
109#pragma omp parallel for
110 for (std::ptrdiff_t j = 0; j < std::ssize(elements); ++j)
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)
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]; });
126 element_node_ids_bulk, connected_element_ids_per_node);
128 if (bulk_element_ids.empty())
130 ERR(
"No element could be found for the subdomain element {:d}. "
131 "Corresponding bulk mesh node ids are:",
133 for (
auto const i : element_node_ids_bulk)
138 "Expect at least one element to be found in the bulk mesh.");
141 bulk_element_ids_map[e->getID()] = std::move(bulk_element_ids);
144 return bulk_element_ids_map;
150 std::vector<std::size_t>
const& values,
154 if (!properties.existsPropertyVector<std::size_t>(property_name,
157 addPropertyToMesh<std::size_t>(mesh, property_name, mesh_item_type, 1,
165 auto& original_property =
167 if (ranges::equal(original_property, values))
170 "There is already a '{:s}' property present in the subdomain mesh "
171 "'{:s}' and it is equal to the newly computed values.",
172 property_name, mesh.
getName());
180 "There is already a '{:s}' property present in the subdomain mesh "
181 "'{:s}' and it is not equal to the newly computed values.",
185 if (!force_overwrite)
187 OGS_FATAL(
"The force overwrite flag was not specified, exiting.");
190 INFO(
"Overwriting '{:s}' property.", property_name);
191 original_property.assign(values);
200 bool const force_overwrite =
false)
204 auto const& bulk_node_ids =
205 identifySubdomainMeshNodes(subdomain_mesh, mesh_node_searcher);
206 INFO(
"identifySubdomainMesh(): identifySubdomainMeshNodes took {:g} s",
209 updateOrCheckExistingSubdomainProperty(
214 auto const& bulk_element_ids =
215 identifySubdomainMeshElements(subdomain_mesh, bulk_mesh);
216 INFO(
"identifySubdomainMesh(): identifySubdomainMeshElements took {:g} s",
224 if (all_of(begin(bulk_element_ids), end(bulk_element_ids),
225 [](std::vector<std::size_t>
const& v) {
return v.size() == 1; }))
229 std::vector<std::size_t> unique_bulk_element_ids;
230 unique_bulk_element_ids.reserve(bulk_element_ids.size());
231 transform(begin(bulk_element_ids), end(bulk_element_ids),
232 back_inserter(unique_bulk_element_ids),
233 [](std::vector<std::size_t>
const& v) {
return v[0]; });
235 updateOrCheckExistingSubdomainProperty(
246 std::vector<std::size_t> flat_bulk_element_ids;
247 flat_bulk_element_ids.reserve(2 * bulk_element_ids.size());
248 std::vector<std::size_t> number_of_bulk_element_ids;
249 number_of_bulk_element_ids.reserve(bulk_element_ids.size());
251 for (std::vector<std::size_t>
const& v : bulk_element_ids)
253 number_of_bulk_element_ids.push_back(v.size());
254 flat_bulk_element_ids.insert(end(flat_bulk_element_ids), begin(v),
258 updateOrCheckExistingSubdomainProperty(
259 subdomain_mesh,
"number_bulk_elements", number_of_bulk_element_ids,
261 updateOrCheckExistingSubdomainProperty(
Definition of the Element class.
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition of the Mesh class.
Definition of the Node class.
Definition of the RunTime class.
double elapsed() const
Get the elapsed time in seconds.
void start()
Start the timer.
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Properties & getProperties()
const std::string getName() const
Get name of the mesh.
std::size_t getNumberOfNodes() const
Get the number of nodes.
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
std::size_t getNumberOfElements() const
Get the number of elements.
PropertyVector< T > const * getPropertyVector(std::string_view name) const
constexpr ranges::views::view_closure ids
For an element of a range view return its id.
constexpr std::string_view getBulkIDString(MeshItemType mesh_item_type)
std::size_t getNodeIndex(Element const &element, unsigned const idx)
PropertyVector< std::size_t > const * bulkNodeIDs(Mesh const &mesh)
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)