Add support for domain of vertices.

src/smith/numerics/functional/domain.cpp

@@ -109,6 +109,41 @@ Domain Domain::ofEdges(const mesh_t& mesh, std::function<bool(std::vector<vec3>)
   return domain_of_edges<3>(mesh, func);
 }
 
+template <int d>
+static Domain domain_of_vertices(const mesh_t& mesh, std::function<bool(tensor<double, d>)> 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<double, d> 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<bool(vec2)> func)
+{
+  return domain_of_vertices<2>(mesh, func);
+}
+
+Domain Domain::ofVertices(const mesh_t& mesh, std::function<bool(vec3)> func)
+{
+  return domain_of_vertices<3>(mesh, func);
+}
+
 ///////////////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////////////
 
@@ -267,7 +302,10 @@ Domain Domain::ofElements(const mesh_t& mesh, std::function<bool(std::vector<vec
 
 void Domain::addElement(int geom_id, int elem_id, mfem::Geometry::Type element_geometry)
 {
-  if (element_geometry == mfem::Geometry::SEGMENT) {
+  if (element_geometry == mfem::Geometry::POINT) {
+    vertex_ids_.push_back(geom_id);
+    mfem_vertex_ids_.push_back(elem_id);
+  } else if (element_geometry == mfem::Geometry::SEGMENT) {
     edge_ids_.push_back(geom_id);
     mfem_edge_ids_.push_back(elem_id);
   } else if (element_geometry == mfem::Geometry::TRIANGLE) {
@@ -478,7 +516,13 @@ mfem::Array<int> Domain::dof_list(const smith::fes_t* fes) const
   }
 
   if (dim_ == 0) {
-    // sam: what to do with vertex sets?
+    GetDofs = [&](int i, mfem::Array<int>& 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);

src/smith/numerics/functional/domain.hpp

@@ -82,12 +82,14 @@ struct Domain {
   ///@}
 
   /// @cond
+  std::vector<int> vertex_ids_;
   std::vector<int> edge_ids_;
   std::vector<int> tri_ids_;
   std::vector<int> quad_ids_;
   std::vector<int> tet_ids_;
   std::vector<int> hex_ids_;
 
+  std::vector<int> mfem_vertex_ids_;
   std::vector<int> mfem_edge_ids_;
   std::vector<int> mfem_tri_ids_;
   std::vector<int> mfem_quad_ids_;
@@ -170,6 +172,7 @@ 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_;
@@ -182,6 +185,7 @@ struct Domain {
   /// @brief get elements by geometry type
   const std::vector<int>& 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<uint32_t, mfem::Geometry::NUM_GEOMETRIES> geometry_counts(cons
 {
   std::array<uint32_t, mfem::Geometry::NUM_GEOMETRIES> counts{};
 
-  constexpr std::array<mfem::Geometry::Type, 5> geometries = {mfem::Geometry::SEGMENT, mfem::Geometry::TRIANGLE,
-                                                              mfem::Geometry::SQUARE, mfem::Geometry::TETRAHEDRON,
-                                                              mfem::Geometry::CUBE};
+  constexpr std::array<mfem::Geometry::Type, 6> 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());
   }

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

@@ -5,6 +5,7 @@
 // SPDX-License-Identifier: (BSD-3-Clause)
 
 #include <algorithm>
+#include <cmath>
 #include <string>
 #include <vector>
 
@@ -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<int> 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;

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<bool(vec2)> 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<bool(vec3)> 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();

src/smith/physics/mesh.hpp

@@ -119,6 +119,14 @@ class Mesh {
   smith::Domain& addDomainOfBodyElements(const std::string& domain_name,
                                          std::function<bool(std::vector<vec2>, 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<bool(vec3)> 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<bool(vec2)> func);
+
   /// @brief get space associated with shape displacement
   const mfem::ParFiniteElementSpace& shapeDisplacementSpace();
 

src/smith/physics/tests/test_mesh_wrapper.cpp

@@ -5,6 +5,7 @@
 // SPDX-License-Identifier: (BSD-3-Clause)
 
 #include <iostream>
+#include <cmath>
 #include <string>
 #include <memory>
 #include <utility>
@@ -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)