From d4d79ebc968d80719bae651dda76cdd220c12152 Mon Sep 17 00:00:00 2001 From: Sam Mish Date: Tue, 26 Nov 2024 15:12:45 -0800 Subject: [PATCH] fix bug when creating domains of boundary elements, delete vertex_ids, address doxygen warnings --- src/serac/numerics/functional/domain.cpp | 46 +------------- src/serac/numerics/functional/domain.hpp | 11 ++-- .../functional/tests/domain_tests.cpp | 60 ------------------- 3 files changed, 8 insertions(+), 109 deletions(-) diff --git a/src/serac/numerics/functional/domain.cpp b/src/serac/numerics/functional/domain.cpp index cfeb5f7b8..55dfe9f81 100644 --- a/src/serac/numerics/functional/domain.cpp +++ b/src/serac/numerics/functional/domain.cpp @@ -48,44 +48,6 @@ std::vector> gather(const mfem::Vector& coordinates, mfem::Arr return x; } -template -static Domain domain_of_vertices(const mfem::Mesh& mesh, std::function)> predicate) -{ - assert(mesh.SpaceDimension() == d); - - Domain output{mesh, 0 /* points are 0-dimensional */}; - - // layout is undocumented, but it seems to be - // [x1, x2, x3, ..., y1, y2, y3 ..., (z1, z2, z3, ...)] - mfem::Vector vertices; - mesh.GetVertices(vertices); - - // vertices that satisfy the predicate are added to the domain - int num_vertices = mesh.GetNV(); - 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.vertex_ids_.push_back(i); - } - } - - return output; -} - -Domain Domain::ofVertices(const mfem::Mesh& mesh, std::function func) -{ - return domain_of_vertices(mesh, func); -} - -Domain Domain::ofVertices(const mfem::Mesh& mesh, std::function func) -{ - return domain_of_vertices(mesh, func); -} - /////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////// @@ -378,13 +340,13 @@ static Domain domain_of_boundary_elems(const mfem::Mesh& break; case mfem::Geometry::TRIANGLE: if (add) { - output.addElement(edge_id, f, geom); + output.addElement(tri_id, f, geom); } tri_id++; break; case mfem::Geometry::SQUARE: if (add) { - output.addElement(edge_id, f, geom); + output.addElement(quad_id, f, geom); } quad_id++; break; @@ -542,10 +504,6 @@ Domain set_operation(set_op op, const Domain& a, const Domain& b) using Ids = std::vector; auto apply_set_op = [&op](const Ids& x, const Ids& y) { return set_operation(op, x, y); }; - if (output.dim_ == 0) { - output.vertex_ids_ = apply_set_op(a.vertex_ids_, b.vertex_ids_); - } - auto fill_output_lists = [apply_set_op, &output](const Ids& a_ids, const Ids& a_mfem_ids, const Ids& b_ids, const Ids& b_mfem_ids, mfem::Geometry::Type g) { auto output_ids = apply_set_op(a_ids, b_ids); diff --git a/src/serac/numerics/functional/domain.hpp b/src/serac/numerics/functional/domain.hpp index 66f6e8ca8..7ecfdcfae 100644 --- a/src/serac/numerics/functional/domain.hpp +++ b/src/serac/numerics/functional/domain.hpp @@ -37,12 +37,10 @@ struct Domain { /// @brief whether the elements in this domain are on the boundary or not Type type_; - std::vector vertex_ids_; - ///@{ /// @name ElementIds /// Indices of elements contained in the domain. - /// The first set, (vertex_ids_, edge_ids_, ...) hold the index of an element in + /// The first set, (edge_ids_, tri_ids, ...) hold the index of an element in /// this Domain in the set of all elements of like geometry in the mesh. /// For example, if edge_ids_[0] = 5, then element 0 in this domain is element /// 5 in the grouping of all edges in the mesh. In other words, these lists @@ -68,6 +66,9 @@ struct Domain { /// methods to add new entities, as this requires you to add both entries and /// keep the corresponding lists in sync. You are discouraged from /// manipulating these lists directly. + ///@} + + /// @cond std::vector edge_ids_; std::vector tri_ids_; std::vector quad_ids_; @@ -79,8 +80,9 @@ struct Domain { std::vector mfem_quad_ids_; std::vector mfem_tet_ids_; std::vector mfem_hex_ids_; - ///@} + /// @endcond + /// @brief construct an "empty" domain, to later be populated later with addElement member functions Domain(const mfem::Mesh& m, int d, Type type = Domain::Type::Elements) : mesh_(m), dim_(d), type_(type) {} /** @@ -143,7 +145,6 @@ 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_; diff --git a/src/serac/numerics/functional/tests/domain_tests.cpp b/src/serac/numerics/functional/tests/domain_tests.cpp index 1fb197dc2..7893107a9 100644 --- a/src/serac/numerics/functional/tests/domain_tests.cpp +++ b/src/serac/numerics/functional/tests/domain_tests.cpp @@ -31,66 +31,6 @@ tensor average(std::vector >& positions) return total / double(positions.size()); } -TEST(domain, of_vertices) -{ - { - auto mesh = import_mesh("onehex.mesh"); - Domain d0 = Domain::ofVertices(mesh, std::function([](vec3 x) { return x[0] < 0.5; })); - EXPECT_EQ(d0.vertex_ids_.size(), 4); - EXPECT_EQ(d0.dim_, 0); - - Domain d1 = Domain::ofVertices(mesh, std::function([](vec3 x) { return x[1] < 0.5; })); - EXPECT_EQ(d1.vertex_ids_.size(), 4); - EXPECT_EQ(d1.dim_, 0); - - Domain d2 = d0 | d1; - EXPECT_EQ(d2.vertex_ids_.size(), 6); - EXPECT_EQ(d2.dim_, 0); - - Domain d3 = d0 & d1; - EXPECT_EQ(d3.vertex_ids_.size(), 2); - EXPECT_EQ(d3.dim_, 0); - } - - { - auto mesh = import_mesh("onetet.mesh"); - Domain d0 = Domain::ofVertices(mesh, std::function([](vec3 x) { return x[0] < 0.5; })); - EXPECT_EQ(d0.vertex_ids_.size(), 3); - EXPECT_EQ(d0.dim_, 0); - - Domain d1 = Domain::ofVertices(mesh, std::function([](vec3 x) { return x[1] < 0.5; })); - EXPECT_EQ(d1.vertex_ids_.size(), 3); - EXPECT_EQ(d1.dim_, 0); - - Domain d2 = d0 | d1; - EXPECT_EQ(d2.vertex_ids_.size(), 4); - EXPECT_EQ(d2.dim_, 0); - - Domain d3 = d0 & d1; - EXPECT_EQ(d3.vertex_ids_.size(), 2); - EXPECT_EQ(d3.dim_, 0); - } - - { - auto mesh = import_mesh("beam-quad.mesh"); - Domain d0 = Domain::ofVertices(mesh, std::function([](vec2 x) { return x[0] < 0.5; })); - EXPECT_EQ(d0.vertex_ids_.size(), 2); - EXPECT_EQ(d0.dim_, 0); - - Domain d1 = Domain::ofVertices(mesh, std::function([](vec2 x) { return x[1] < 0.5; })); - EXPECT_EQ(d1.vertex_ids_.size(), 9); - EXPECT_EQ(d1.dim_, 0); - - Domain d2 = d0 | d1; - EXPECT_EQ(d2.vertex_ids_.size(), 10); - EXPECT_EQ(d2.dim_, 0); - - Domain d3 = d0 & d1; - EXPECT_EQ(d3.vertex_ids_.size(), 1); - EXPECT_EQ(d3.dim_, 0); - } -} - TEST(domain, of_edges) { {