From 4f0e2b9d616b630820c2000b206ab2d6eef9b649 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Mon, 23 Mar 2026 14:06:16 -0700 Subject: [PATCH] Add support for domain of vertices. --- src/smith/numerics/functional/domain.cpp | 52 ++++++++++++++++++- src/smith/numerics/functional/domain.hpp | 15 ++++-- .../functional/tests/domain_tests.cpp | 27 ++++++++++ src/smith/physics/mesh.cpp | 14 +++++ src/smith/physics/mesh.hpp | 8 +++ src/smith/physics/tests/test_mesh_wrapper.cpp | 5 ++ 6 files changed, 114 insertions(+), 7 deletions(-) diff --git a/src/smith/numerics/functional/domain.cpp b/src/smith/numerics/functional/domain.cpp index 9acbe708c4..2cce995843 100644 --- a/src/smith/numerics/functional/domain.cpp +++ b/src/smith/numerics/functional/domain.cpp @@ -109,6 +109,41 @@ Domain Domain::ofEdges(const mesh_t& mesh, std::function) return domain_of_edges<3>(mesh, func); } +template +static Domain domain_of_vertices(const mesh_t& mesh, std::function)> predicate) +{ + assert(mesh.SpaceDimension() == d); + + Domain output{mesh, 0, Domain::Type::Elements}; + + mfem::Vector vertices; + mesh.GetVertices(vertices); + + int num_vertices = vertices.Size() / d; + for (int i = 0; i < num_vertices; i++) { + tensor x{}; + for (int j = 0; j < d; j++) { + x[j] = vertices[j * num_vertices + i]; + } + + if (predicate(x)) { + output.addElement(i, i, mfem::Geometry::POINT); + } + } + + return output; +} + +Domain Domain::ofVertices(const mesh_t& mesh, std::function func) +{ + return domain_of_vertices<2>(mesh, func); +} + +Domain Domain::ofVertices(const mesh_t& mesh, std::function func) +{ + return domain_of_vertices<3>(mesh, func); +} + /////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////// @@ -267,7 +302,10 @@ Domain Domain::ofElements(const mesh_t& mesh, std::function Domain::dof_list(const smith::fes_t* fes) const } if (dim_ == 0) { - // sam: what to do with vertex sets? + GetDofs = [&](int i, mfem::Array& vdofs) { return fes->GetVertexDofs(i, vdofs); }; + for (auto vertex_id : mfem_vertex_ids_) { + GetDofs(vertex_id, elem_dofs); + for (int i = 0; i < elem_dofs.Size(); i++) { + dof_ids.insert(elem_dofs[i]); + } + } } if (dim_ == 1) { @@ -743,6 +787,10 @@ Domain set_operation(SET_OPERATION op, const Domain& a, const Domain& b) fill_combined_lists(a.edge_ids_, a.mfem_edge_ids_, b.edge_ids_, b.mfem_edge_ids_, mfem::Geometry::SEGMENT); } + if (combined.dim_ == 0) { + fill_combined_lists(a.vertex_ids_, a.mfem_vertex_ids_, b.vertex_ids_, b.mfem_vertex_ids_, mfem::Geometry::POINT); + } + if (combined.dim_ == 2) { fill_combined_lists(a.tri_ids_, a.mfem_tri_ids_, b.tri_ids_, b.mfem_tri_ids_, mfem::Geometry::TRIANGLE); fill_combined_lists(a.quad_ids_, a.mfem_quad_ids_, b.quad_ids_, b.mfem_quad_ids_, mfem::Geometry::SQUARE); diff --git a/src/smith/numerics/functional/domain.hpp b/src/smith/numerics/functional/domain.hpp index a7fe5747e0..791849497b 100644 --- a/src/smith/numerics/functional/domain.hpp +++ b/src/smith/numerics/functional/domain.hpp @@ -82,12 +82,14 @@ struct Domain { ///@} /// @cond + std::vector vertex_ids_; std::vector edge_ids_; std::vector tri_ids_; std::vector quad_ids_; std::vector tet_ids_; std::vector hex_ids_; + std::vector mfem_vertex_ids_; std::vector mfem_edge_ids_; std::vector mfem_tri_ids_; std::vector mfem_quad_ids_; @@ -170,6 +172,7 @@ struct Domain { /// @brief get elements by geometry type const std::vector& get(mfem::Geometry::Type geom) const { + if (geom == mfem::Geometry::POINT) return vertex_ids_; if (geom == mfem::Geometry::SEGMENT) return edge_ids_; if (geom == mfem::Geometry::TRIANGLE) return tri_ids_; if (geom == mfem::Geometry::SQUARE) return quad_ids_; @@ -182,6 +185,7 @@ struct Domain { /// @brief get elements by geometry type const std::vector& get_mfem_ids(mfem::Geometry::Type geom) const { + if (geom == mfem::Geometry::POINT) return mfem_vertex_ids_; if (geom == mfem::Geometry::SEGMENT) return mfem_edge_ids_; if (geom == mfem::Geometry::TRIANGLE) return mfem_tri_ids_; if (geom == mfem::Geometry::SQUARE) return mfem_quad_ids_; @@ -196,7 +200,8 @@ struct Domain { */ int total_elements() const { - return int(edge_ids_.size() + tri_ids_.size() + quad_ids_.size() + tet_ids_.size() + hex_ids_.size()); + return int(vertex_ids_.size() + edge_ids_.size() + tri_ids_.size() + quad_ids_.size() + tet_ids_.size() + + hex_ids_.size()); } /** @@ -209,7 +214,7 @@ struct Domain { int total = 0; offsets[mfem::Geometry::POINT] = total; - total += 0; // vertices; + total += int(vertex_ids_.size()); offsets[mfem::Geometry::SEGMENT] = total; total += int(edge_ids_.size()); offsets[mfem::Geometry::TRIANGLE] = total; @@ -301,9 +306,9 @@ inline std::array geometry_counts(cons { std::array counts{}; - constexpr std::array geometries = {mfem::Geometry::SEGMENT, mfem::Geometry::TRIANGLE, - mfem::Geometry::SQUARE, mfem::Geometry::TETRAHEDRON, - mfem::Geometry::CUBE}; + constexpr std::array geometries = {mfem::Geometry::POINT, mfem::Geometry::SEGMENT, + mfem::Geometry::TRIANGLE, mfem::Geometry::SQUARE, + mfem::Geometry::TETRAHEDRON, mfem::Geometry::CUBE}; for (auto geom : geometries) { counts[uint32_t(geom)] = uint32_t(domain.get(geom).size()); } diff --git a/src/smith/numerics/functional/tests/domain_tests.cpp b/src/smith/numerics/functional/tests/domain_tests.cpp index a73da97f0e..ad9bf51f03 100644 --- a/src/smith/numerics/functional/tests/domain_tests.cpp +++ b/src/smith/numerics/functional/tests/domain_tests.cpp @@ -5,6 +5,7 @@ // SPDX-License-Identifier: (BSD-3-Clause) #include +#include #include #include @@ -290,6 +291,32 @@ TEST(domain, of_elements) } } +TEST(domain, of_vertices_finds_dofs) +{ + constexpr int dim = 2; + constexpr int p = 2; + auto bmesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto mesh = smith::mesh::refineAndDistribute(std::move(bmesh)); + + auto fec = mfem::H1_FECollection(p, dim); + auto fes = mfem::ParFiniteElementSpace(mesh.get(), &fec); + + Domain lower_left = + Domain::ofVertices(*mesh, [](vec2 x) { return std::abs(x[0]) < 1e-12 && std::abs(x[1]) < 1e-12; }); + EXPECT_EQ(lower_left.dim_, 0); + EXPECT_EQ(lower_left.vertex_ids_.size(), 1); + + mfem::Array dof_indices = lower_left.dof_list(&fes); + EXPECT_EQ(dof_indices.Size(), 1); + + Domain upper_right = + Domain::ofVertices(*mesh, [](vec2 x) { return std::abs(x[0] - 1.0) < 1e-12 && std::abs(x[1] - 1.0) < 1e-12; }); + EXPECT_EQ(upper_right.vertex_ids_.size(), 1); + + dof_indices = (lower_left | upper_right).dof_list(&fes); + EXPECT_EQ(dof_indices.Size(), 2); +} + TEST(domain, entireDomain2d) { constexpr int dim = 2; diff --git a/src/smith/physics/mesh.cpp b/src/smith/physics/mesh.cpp index 83ed0b5257..6be7b463d0 100644 --- a/src/smith/physics/mesh.cpp +++ b/src/smith/physics/mesh.cpp @@ -143,6 +143,20 @@ smith::Domain& Mesh::addDomainOfBodyElements(const std::string& domain_name, return domain(domain_name); } +smith::Domain& Mesh::addDomainOfVertices(const std::string& domain_name, std::function func) +{ + errorIfDomainExists(domain_name); + domains_.emplace(domain_name, Domain::ofVertices(*mfem_mesh_, func)); + return domain(domain_name); +} + +smith::Domain& Mesh::addDomainOfVertices(const std::string& domain_name, std::function func) +{ + errorIfDomainExists(domain_name); + domains_.emplace(domain_name, Domain::ofVertices(*mfem_mesh_, func)); + return domain(domain_name); +} + const mfem::ParFiniteElementSpace& Mesh::shapeDisplacementSpace() { return smith::StateManager::shapeDisplacement(tag()).space(); diff --git a/src/smith/physics/mesh.hpp b/src/smith/physics/mesh.hpp index cac8efe4ce..a32794fc50 100644 --- a/src/smith/physics/mesh.hpp +++ b/src/smith/physics/mesh.hpp @@ -119,6 +119,14 @@ class Mesh { smith::Domain& addDomainOfBodyElements(const std::string& domain_name, std::function, int)> func); + /// @brief create domain of vertices with specified name on a 3D mesh + /// The second argument is a function taking the nodal coordinate of the vertex + smith::Domain& addDomainOfVertices(const std::string& domain_name, std::function func); + + /// @brief create domain of vertices with specified name on a 2D mesh + /// The second argument is a function taking the nodal coordinate of the vertex + smith::Domain& addDomainOfVertices(const std::string& domain_name, std::function func); + /// @brief get space associated with shape displacement const mfem::ParFiniteElementSpace& shapeDisplacementSpace(); diff --git a/src/smith/physics/tests/test_mesh_wrapper.cpp b/src/smith/physics/tests/test_mesh_wrapper.cpp index 913d3bf506..45ede4e939 100644 --- a/src/smith/physics/tests/test_mesh_wrapper.cpp +++ b/src/smith/physics/tests/test_mesh_wrapper.cpp @@ -5,6 +5,7 @@ // SPDX-License-Identifier: (BSD-3-Clause) #include +#include #include #include #include @@ -58,6 +59,10 @@ TEST(Mesh, TestMeshFromSerialMesh) }; mesh->addDomainOfBodyElements("left", is_in_left); EXPECT_EQ(mesh->domain("left").total_elements(), 4); + + mesh->addDomainOfVertices( + "origin", [](vec3 x) { return std::abs(x[0]) < 1e-12 && std::abs(x[1]) < 1e-12 && std::abs(x[2]) < 1e-12; }); + EXPECT_EQ(mesh->domain("origin").total_elements(), 1); } TEST(Mesh, TestMeshFromFile)