Merge pull request #1286 from LLNL/talamini/bugfix/domain_element_type

Fix incorrect type in some boolean operations on Domain

src/serac/numerics/functional/domain.cpp

@@ -56,7 +56,7 @@ static Domain domain_of_edges(const mesh_t& mesh, std::function<T> predicate)
 {
   assert(mesh.SpaceDimension() == d);
 
-  Domain output{mesh, 1 /* edges are 1-dimensional */};
+  Domain output{mesh, 1 /* edges are 1-dimensional */, Domain::Type::Elements};
 
   // layout is undocumented, but it seems to be
   // [x1, x2, x3, ..., y1, y2, y3 ..., (z1, z2, z3, ...)]
@@ -109,7 +109,7 @@ static Domain domain_of_faces(const mesh_t& mesh, std::function<bool(std::vector
 {
   assert(mesh.SpaceDimension() == d);
 
-  Domain output{mesh, 2 /* faces are 2-dimensional */};
+  Domain output{mesh, 2 /* faces are 2-dimensional */, Domain::Type::Elements};
 
   // layout is undocumented, but it seems to be
   // [x1, x2, x3, ..., y1, y2, y3 ..., (z1, z2, z3, ...)]
@@ -189,7 +189,7 @@ static Domain domain_of_elems(const mesh_t& mesh, std::function<bool(std::vector
 {
   assert(mesh.SpaceDimension() == d);
 
-  Domain output{mesh, mesh.SpaceDimension() /* elems can be 2 or 3 dimensional */};
+  Domain output{mesh, mesh.SpaceDimension() /* elems can be 2 or 3 dimensional */, Domain::Type::Elements};
 
   // layout is undocumented, but it seems to be
   // [x1, x2, x3, ..., y1, y2, y3 ..., (z1, z2, z3, ...)]
@@ -582,8 +582,9 @@ Domain set_operation(SET_OPERATION op, const Domain& a, const Domain& b)
 {
   assert(&a.mesh_ == &b.mesh_);
   assert(a.dim_ == b.dim_);
+  assert(a.type_ == b.type_);
 
-  Domain combined{a.mesh_, a.dim_};
+  Domain combined{a.mesh_, a.dim_, a.type_};
 
   using Ids         = std::vector<int>;
   auto apply_set_op = [&op](const Ids& x, const Ids& y) { return set_operation(op, x, y); };

src/serac/numerics/functional/domain.hpp

@@ -99,7 +99,7 @@ struct Domain {
   /**
    * @brief empty Domain constructor, with connectivity info to be populated later
    */
-  Domain(const mesh_t& m, int d, Type type = Domain::Type::Elements) : mesh_(m), dim_(d), type_(type) {}
+  Domain(const mesh_t& m, int d, Type type) : mesh_(m), dim_(d), type_(type) {}
 
   /**
    * @brief create a domain from some subset of the vertices in an mfem::Mesh

src/serac/numerics/functional/tests/domain_tests.cpp

@@ -394,6 +394,74 @@ TEST(domain, of3dElementsFindsDofs)
   EXPECT_EQ(dof_indices.Size(), 113);
 }
 
+TEST(domain, of2dBoundaryElementsFindsDofs)
+{
+  constexpr int dim  = 2;
+  constexpr int p    = 2;
+  auto          mesh = import_mesh("patch2D_tris_and_quads.mesh");
+
+  auto find_right_boundary = [](std::vector<vec2> vertices, int /* attr */) {
+    return std::all_of(vertices.begin(), vertices.end(), [](vec2 X) { return X[0] > 1.0 - 1e-2; });
+  };
+
+  Domain d0 = Domain::ofBoundaryElements(mesh, find_right_boundary);
+  EXPECT_EQ(d0.edge_ids_.size(), 1);
+
+  auto fec = mfem::H1_FECollection(p, dim);
+  auto fes = mfem::FiniteElementSpace(&mesh, &fec);
+
+  mfem::Array<int> dof_indices = d0.dof_list(&fes);
+
+  EXPECT_EQ(dof_indices.Size(), 3);
+
+  auto find_top_boundary = [](std::vector<vec2> vertices, int /* attr */) {
+    return std::all_of(vertices.begin(), vertices.end(), [](vec2 X) { return X[1] > 1.0 - 1e-2; });
+  };
+
+  Domain d1 = Domain::ofBoundaryElements(mesh, find_top_boundary);
+  EXPECT_EQ(d1.edge_ids_.size(), 1);
+
+  Domain d2 = d0 | d1;
+
+  dof_indices = d2.dof_list(&fes);
+
+  EXPECT_EQ(dof_indices.Size(), 5);
+}
+
+TEST(domain, of3dBoundaryElementsFindsDofs)
+{
+  constexpr int dim  = 3;
+  constexpr int p    = 2;
+  auto          mesh = import_mesh("patch3D_tets.mesh");
+
+  auto find_xmax_boundary = [](std::vector<vec3> vertices, int /* attr */) {
+    return std::all_of(vertices.begin(), vertices.end(), [](vec3 X) { return X[0] > 1.0 - 1e-2; });
+  };
+
+  Domain d0 = Domain::ofBoundaryElements(mesh, find_xmax_boundary);
+  EXPECT_EQ(d0.tri_ids_.size(), 2);
+
+  auto fec = mfem::H1_FECollection(p, dim);
+  auto fes = mfem::FiniteElementSpace(&mesh, &fec);
+
+  mfem::Array<int> dof_indices = d0.dof_list(&fes);
+
+  EXPECT_EQ(dof_indices.Size(), 9);
+
+  auto find_ymax_boundary = [](std::vector<vec3> vertices, int /* attr */) {
+    return std::all_of(vertices.begin(), vertices.end(), [](vec3 X) { return X[1] > 1.0 - 1e-2; });
+  };
+
+  Domain d1 = Domain::ofBoundaryElements(mesh, find_ymax_boundary);
+  EXPECT_EQ(d1.tri_ids_.size(), 2);
+
+  Domain d2 = d0 | d1;
+
+  dof_indices = d2.dof_list(&fes);
+
+  EXPECT_EQ(dof_indices.Size(), 15);
+}
+
 int main(int argc, char* argv[])
 {
   int num_procs, myid;