Skip to content

Commit

Permalink
fix bug when creating domains of boundary elements, delete vertex_ids…
Browse files Browse the repository at this point in the history 
…, address doxygen warnings
  • Loading branch information
@samuelpmishLLNL
samuelpmishLLNL committed Nov 26, 2024
1 parent 9cb9051 commit d4d79eb
Show file tree
Hide file tree
Showing 3 changed files with 8 additions and 109 deletions.
46 changes: 2 additions & 44 deletions src/serac/numerics/functional/domain.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,44 +48,6 @@ std::vector<tensor<double, d>> gather(const mfem::Vector& coordinates, mfem::Arr
return x;
}

template <int d>
static Domain domain_of_vertices(const mfem::Mesh& mesh, std::function<bool(tensor<double, d>)> 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<double, d> 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<bool(vec2)> func)
{
return domain_of_vertices(mesh, func);
}

Domain Domain::ofVertices(const mfem::Mesh& mesh, std::function<bool(vec3)> func)
{
return domain_of_vertices(mesh, func);
}

///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////

Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -542,10 +504,6 @@ Domain set_operation(set_op op, const Domain& a, const Domain& b)
using Ids = std::vector<int>;
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);
Expand Down
11 changes: 6 additions & 5 deletions src/serac/numerics/functional/domain.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,12 +37,10 @@ struct Domain {
/// @brief whether the elements in this domain are on the boundary or not
Type type_;

std::vector<int> 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
Expand All @@ -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<int> edge_ids_;
std::vector<int> tri_ids_;
std::vector<int> quad_ids_;
Expand All @@ -79,8 +80,9 @@ struct Domain {
std::vector<int> mfem_quad_ids_;
std::vector<int> mfem_tet_ids_;
std::vector<int> 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) {}

/**
Expand Down Expand Up @@ -143,7 +145,6 @@ struct Domain {
/// @brief get elements by geometry type
const std::vector<int>& 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_;
Expand Down
60 changes: 0 additions & 60 deletions src/serac/numerics/functional/tests/domain_tests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,66 +31,6 @@ tensor<double, dim> average(std::vector<tensor<double, dim> >& 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)
{
{
Expand Down

0 comments on commit d4d79eb

Please sign in to comment.