diff --git a/Applications/Utils/MeshGeoTools/IdentifySubdomains.cpp b/Applications/Utils/MeshGeoTools/IdentifySubdomains.cpp
index 8d1d3e592945112e8f7f8c6b0afe77c1fdba4b66..c97e09e7a982b09214128e0f34b63dbcb3f42a47 100644
--- a/Applications/Utils/MeshGeoTools/IdentifySubdomains.cpp
+++ b/Applications/Utils/MeshGeoTools/IdentifySubdomains.cpp
@@ -26,7 +26,12 @@ std::vector<std::unique_ptr<MeshLib::Mesh>> readMeshes(
for (auto const& filename : filenames)
{
- meshes.emplace_back(MeshLib::IO::readMeshFromFile(filename));
+ auto mesh = MeshLib::IO::readMeshFromFile(filename);
+ if (mesh == nullptr)
+ {
+ OGS_FATAL("Could not read mesh from '%s' file.", filename.c_str());
+ }
+ meshes.emplace_back(mesh);
}
if (meshes.empty())
{
@@ -95,6 +100,11 @@ int main(int argc, char* argv[])
//
std::unique_ptr<MeshLib::Mesh> bulk_mesh{
MeshLib::IO::readMeshFromFile(bulk_mesh_arg.getValue())};
+ if (bulk_mesh == nullptr)
+ {
+ OGS_FATAL("Could not read bulk mesh from '%s'",
+ bulk_mesh_arg.getValue().c_str());
+ }
//
// Read the subdomain meshes.
diff --git a/Applications/Utils/Tests.cmake b/Applications/Utils/Tests.cmake
index d99d4f39690c4f2cfb75129041bbc9daf5e68e6e..bb69b4cbcc44e6baaf46ecdd03831c661eda467c 100644
--- a/Applications/Utils/Tests.cmake
+++ b/Applications/Utils/Tests.cmake
@@ -68,6 +68,32 @@ AddTest(
2D_mesh_bottom.vtu check_2D_mesh_bottom.vtu bulk_element_ids bulk_element_ids 0 0
)
+AddTest(
+ NAME identifySubdomains_riverTriangleMesh
+ PATH MeshGeoToolsLib/IdentifySubdomains
+ EXECUTABLE identifySubdomains
+ EXECUTABLE_ARGS -m river_domain_triangle.vtu -o ${Data_BINARY_DIR}/MeshGeoToolsLib/IdentifySubdomains/triangle_ -- river_bc.vtu
+ REQUIREMENTS NOT OGS_USE_MPI
+ TESTER vtkdiff
+ DIFF_DATA
+ river_bc_triangle.vtu triangle_river_bc.vtu bulk_node_ids bulk_node_ids 0 0
+ #river_bc_triangle.vtu triangle_river_bc.vtu bulk_element_ids bulk_element_ids 0 0 # TODO (naumov) Needs extension of vtkdiff to FieldData
+ river_bc_triangle.vtu triangle_river_bc.vtu number_bulk_elements number_bulk_elements 0 0
+)
+
+AddTest(
+ NAME identifySubdomains_riverPrismMesh
+ PATH MeshGeoToolsLib/IdentifySubdomains
+ EXECUTABLE identifySubdomains
+ EXECUTABLE_ARGS -s 1e-3 -m river_domain_prism.vtu -o ${Data_BINARY_DIR}/MeshGeoToolsLib/IdentifySubdomains/prism_ -- river_bc.vtu
+ REQUIREMENTS NOT OGS_USE_MPI
+ TESTER vtkdiff
+ DIFF_DATA
+ river_bc_prism.vtu prism_river_bc.vtu bulk_node_ids bulk_node_ids 0 0
+ #river_bc_prism.vtu prism_river_bc.vtu bulk_element_ids bulk_element_ids 0 0 # TODO (naumov) Needs extension of vtkdiff to FieldData
+ river_bc_prism.vtu prism_river_bc.vtu number_bulk_elements number_bulk_elements 0 0
+)
+
# Mac is producing slightly different partitioning, so the results are not
# comparable.
AddTest(
diff --git a/MeshGeoToolsLib/IdentifySubdomainMesh.cpp b/MeshGeoToolsLib/IdentifySubdomainMesh.cpp
index e69a202d4b8ca26db3fe3fb64bb42c4ca287641b..5b69d073d90ad38f28af3c494ff5ccdda75bb40b 100644
--- a/MeshGeoToolsLib/IdentifySubdomainMesh.cpp
+++ b/MeshGeoToolsLib/IdentifySubdomainMesh.cpp
@@ -92,19 +92,19 @@ std::vector<std::size_t> findElementsInMesh(
return element_ids;
}
-/// Tries to find unique elements' ids in the bulk mesh.
+/// Tries to find all elements' ids in the bulk mesh.
/// The subdomain elements' nodes are mapped to the bulk mesh and are then used
/// to identify the bulk mesh elements.
-std::vector<std::size_t> identifySubdomainMeshElements(
+std::vector<std::vector<std::size_t>> identifySubdomainMeshElements(
MeshLib::Mesh const& subdomain_mesh, MeshLib::Mesh const& bulk_mesh)
{
auto& properties = subdomain_mesh.getProperties();
auto const& bulk_node_ids =
*properties.getPropertyVector<std::size_t>("bulk_node_ids");
- // Initialize with very large value
- std::vector<std::size_t> bulk_element_ids_map(
- subdomain_mesh.getNumberOfElements(), -1);
+ // Allocate space for all elements for random insertion.
+ std::vector<std::vector<std::size_t>> bulk_element_ids_map(
+ subdomain_mesh.getNumberOfElements());
for (auto* const e : subdomain_mesh.getElements())
{
@@ -120,7 +120,7 @@ std::vector<std::size_t> identifySubdomainMeshElements(
[&bulk_node_ids](std::size_t const id) {
return bulk_node_ids[id];
});
- auto const& bulk_element_ids =
+ std::vector<std::size_t> bulk_element_ids =
findElementsInMesh(bulk_mesh, element_node_ids_bulk);
if (bulk_element_ids.empty())
@@ -131,23 +131,10 @@ std::vector<std::size_t> identifySubdomainMeshElements(
for (auto const i : element_node_ids_bulk)
ERR("\t%d", i);
OGS_FATAL(
- "Expect exactly one element to be found in the bulk mesh.");
- }
- if (bulk_element_ids.size() != 1)
- {
- ERR("%d elements found for the subdomain element %d. "
- "Corresponding bulk mesh node ids are:",
- bulk_element_ids.size(), e->getID());
- for (auto const i : element_node_ids_bulk)
- ERR("\t%d", i);
- ERR("The found bulk mesh element ids are:");
- for (auto const i : bulk_element_ids)
- ERR("\t%d", i);
- OGS_FATAL(
- "Expect exactly one element to be found in the bulk mesh.");
+ "Expect at least one element to be found in the bulk mesh.");
}
- bulk_element_ids_map[e->getID()] = bulk_element_ids.front();
+ bulk_element_ids_map[e->getID()] = std::move(bulk_element_ids);
}
return bulk_element_ids_map;
@@ -156,8 +143,8 @@ std::vector<std::size_t> identifySubdomainMeshElements(
/// Updates or checks the existing mesh's property with the given values.
void updateOrCheckExistingSubdomainProperty(
MeshLib::Mesh& mesh, std::string const& property_name,
- std::vector<std::size_t> values, MeshLib::MeshItemType const mesh_item_type,
- bool const force_overwrite)
+ std::vector<std::size_t> const& values,
+ MeshLib::MeshItemType const mesh_item_type, bool const force_overwrite)
{
auto& properties = mesh.getProperties();
if (!properties.existsPropertyVector<std::size_t>(property_name))
@@ -218,8 +205,49 @@ void identifySubdomainMesh(MeshLib::Mesh& subdomain_mesh,
auto const& bulk_element_ids =
identifySubdomainMeshElements(subdomain_mesh, bulk_mesh);
- updateOrCheckExistingSubdomainProperty(
- subdomain_mesh, "bulk_element_ids", bulk_element_ids,
- MeshLib::MeshItemType::Cell, force_overwrite);
+ // The bulk_element_ids could be of two types: one element per entry---this
+ // is the expected case for the boundary meshes; multiple elements per
+ // entry---this happens if the subdomain mesh lies inside the bulk mesh and
+ // has lower dimension.
+ // First find out the type, then add/check the CellData or FieldData.
+ if (all_of(begin(bulk_element_ids), end(bulk_element_ids),
+ [](std::vector<std::size_t> const& v) { return v.size() == 1; }))
+ {
+ // All vectors are of size 1, so the data can be flattened and
+ // stored in CellData or compared to existing CellData.
+ std::vector<std::size_t> unique_bulk_element_ids;
+ unique_bulk_element_ids.reserve(bulk_element_ids.size());
+ transform(begin(bulk_element_ids), end(bulk_element_ids),
+ back_inserter(unique_bulk_element_ids),
+ [](std::vector<std::size_t> const& v) { return v[0]; });
+
+ updateOrCheckExistingSubdomainProperty(
+ subdomain_mesh, "bulk_element_ids", unique_bulk_element_ids,
+ MeshLib::MeshItemType::Cell, force_overwrite);
+ }
+ else
+ {
+ // Some of the boundary elements are connected to multiple bulk
+ // elements; Store the array in FieldData with additional CellData array
+ // for the number of elements, which also provides the offsets.
+ std::vector<std::size_t> flat_bulk_element_ids;
+ flat_bulk_element_ids.reserve(2 * bulk_element_ids.size()); // Guess.
+ std::vector<std::size_t> number_of_bulk_element_ids;
+ number_of_bulk_element_ids.reserve(bulk_element_ids.size());
+
+ for (std::vector<std::size_t> const& v : bulk_element_ids)
+ {
+ number_of_bulk_element_ids.push_back(v.size());
+ flat_bulk_element_ids.insert(end(flat_bulk_element_ids), begin(v),
+ end(v));
+ }
+
+ updateOrCheckExistingSubdomainProperty(
+ subdomain_mesh, "number_bulk_elements", number_of_bulk_element_ids,
+ MeshLib::MeshItemType::Cell, force_overwrite);
+ updateOrCheckExistingSubdomainProperty(
+ subdomain_mesh, "bulk_element_ids", flat_bulk_element_ids,
+ MeshLib::MeshItemType::IntegrationPoint, force_overwrite);
+ }
}
} // namespace MeshGeoToolsLib
diff --git a/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/river_bc.vtu b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/river_bc.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..2fed6bda560eb2a4cd4687ee94587d445eda7377
--- /dev/null
+++ b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/river_bc.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:61f0e4026ac959f245641a37750b2841e939b58e2aaf1f8f3f755367f0f6538b
+size 6468
diff --git a/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/river_bc_prism.vtu b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/river_bc_prism.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..76457f6eb3ca1aed803e646a3dcc842b800b8641
--- /dev/null
+++ b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/river_bc_prism.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:e417a0e362b6bc610ca502cd2963732ced8cc9efb0b59547a722cbf609d3d528
+size 10012
diff --git a/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/river_bc_triangle.vtu b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/river_bc_triangle.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..dfbe743f9bfae9bd4073246ef0ae1aeea354ad6f
--- /dev/null
+++ b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/river_bc_triangle.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:71248e74e69d1ebc99848fd5a61141e8ecf1155da3d5dbcab291f13b8b15f40b
+size 10012
diff --git a/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/river_domain_prism.vtu b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/river_domain_prism.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..08d436bfe86783cc185536786a6b599bf82b949d
--- /dev/null
+++ b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/river_domain_prism.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:0231b95acb47bf2a99ca4f80a35bbe7779a236e8d20a18c1078f72178e39afaa
+size 61970
diff --git a/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/river_domain_triangle.vtu b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/river_domain_triangle.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..55033e83d26c4b7b3b02da14430eebb0059de97a
--- /dev/null
+++ b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/river_domain_triangle.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5d0fcd0674b781995afa3837ee6156ebd43924fa6c9627ce358b3c9ceb414280
+size 18616