diff --git a/src/serac/numerics/functional/domain.cpp b/src/serac/numerics/functional/domain.cpp index 3d39e2d0a..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); -} - /////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////// @@ -117,13 +79,11 @@ static Domain domain_of_edges(const mfem::Mesh& mesh, std::function predicate int bdr_id = edge_id_to_bdr_id[i]; int attr = (bdr_id > 0) ? mesh.GetBdrAttribute(bdr_id) : -1; if (predicate(x, attr)) { - output.edge_ids_.push_back(i); - output.mfem_edge_ids_.push_back(i); + output.addElement(i, i, mfem::Geometry::SEGMENT); } } else { if (predicate(x)) { - output.edge_ids_.push_back(i); - output.mfem_edge_ids_.push_back(i); + output.addElement(i, i, mfem::Geometry::SEGMENT); } } } @@ -194,12 +154,10 @@ static Domain domain_of_faces(const mfem::Mesh& if (predicate(x, attr)) { if (x.size() == 3) { - output.tri_ids_.push_back(tri_id); - output.mfem_tri_ids_.push_back(i); + output.addElement(tri_id, i, mfem::Geometry::TRIANGLE); } if (x.size() == 4) { - output.quad_ids_.push_back(quad_id); - output.mfem_quad_ids_.push_back(i); + output.addElement(quad_id, i, mfem::Geometry::SQUARE); } } @@ -258,31 +216,27 @@ static Domain domain_of_elems(const mfem::Mesh& switch (x.size()) { case 3: if (add) { - output.tri_ids_.push_back(tri_id); - output.mfem_tri_ids_.push_back(i); + output.addElement(tri_id, i, mfem::Geometry::TRIANGLE); } tri_id++; break; case 4: if constexpr (d == 2) { if (add) { - output.quad_ids_.push_back(quad_id); - output.mfem_quad_ids_.push_back(i); + output.addElement(quad_id, i, mfem::Geometry::SQUARE); } quad_id++; } if constexpr (d == 3) { if (add) { - output.tet_ids_.push_back(tet_id); - output.mfem_tet_ids_.push_back(i); + output.addElement(tet_id, i, mfem::Geometry::TETRAHEDRON); } tet_id++; } break; case 8: if (add) { - output.hex_ids_.push_back(hex_id); - output.mfem_hex_ids_.push_back(i); + output.addElement(hex_id, i, mfem::Geometry::CUBE); } hex_id++; break; @@ -305,6 +259,39 @@ Domain Domain::ofElements(const mfem::Mesh& mesh, std::function(mesh, func); } +void Domain::addElement(int geom_id, int elem_id, mfem::Geometry::Type element_geometry) +{ + 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) { + tri_ids_.push_back(geom_id); + mfem_tri_ids_.push_back(elem_id); + } else if (element_geometry == mfem::Geometry::SQUARE) { + quad_ids_.push_back(geom_id); + mfem_quad_ids_.push_back(elem_id); + } else if (element_geometry == mfem::Geometry::TETRAHEDRON) { + tet_ids_.push_back(geom_id); + mfem_tet_ids_.push_back(elem_id); + } else if (element_geometry == mfem::Geometry::CUBE) { + hex_ids_.push_back(geom_id); + mfem_hex_ids_.push_back(elem_id); + } else { + SLIC_ERROR("unsupported element type"); + } +} + +void Domain::addElements(const std::vector& geom_ids, const std::vector& elem_ids, + mfem::Geometry::Type element_geometry) +{ + SLIC_ERROR_IF(geom_ids.size() != elem_ids.size(), + "To add elements, you must specify a geom_id AND an elem_id for each element"); + + for (std::vector::size_type i = 0; i < geom_ids.size(); ++i) { + addElement(geom_ids[i], elem_ids[i], element_geometry); + } +} + /////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////// @@ -347,22 +334,19 @@ static Domain domain_of_boundary_elems(const mfem::Mesh& switch (geom) { case mfem::Geometry::SEGMENT: if (add) { - output.edge_ids_.push_back(edge_id); - output.mfem_edge_ids_.push_back(f); + output.addElement(edge_id, f, geom); } edge_id++; break; case mfem::Geometry::TRIANGLE: if (add) { - output.tri_ids_.push_back(tri_id); - output.mfem_tri_ids_.push_back(f); + output.addElement(tri_id, f, geom); } tri_id++; break; case mfem::Geometry::SQUARE: if (add) { - output.quad_ids_.push_back(quad_id); - output.mfem_quad_ids_.push_back(f); + output.addElement(quad_id, f, geom); } quad_id++; break; @@ -385,6 +369,9 @@ Domain Domain::ofBoundaryElements(const mfem::Mesh& mesh, std::function(mesh, func); } +/////////////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////////// + mfem::Array Domain::dof_list(mfem::FiniteElementSpace* fes) const { std::set dof_ids; @@ -460,71 +447,32 @@ mfem::Array Domain::dof_list(mfem::FiniteElementSpace* fes) const Domain EntireDomain(const mfem::Mesh& mesh) { - Domain output{mesh, mesh.SpaceDimension() /* elems can be 2 or 3 dimensional */}; - - int tri_id = 0; - int quad_id = 0; - int tet_id = 0; - int hex_id = 0; - - // faces that satisfy the predicate are added to the domain - int num_elems = mesh.GetNE(); - for (int i = 0; i < num_elems; i++) { - auto geom = mesh.GetElementGeometry(i); - - switch (geom) { - case mfem::Geometry::TRIANGLE: - output.tri_ids_.push_back(tri_id++); - break; - case mfem::Geometry::SQUARE: - output.quad_ids_.push_back(quad_id++); - break; - case mfem::Geometry::TETRAHEDRON: - output.tet_ids_.push_back(tet_id++); - break; - case mfem::Geometry::CUBE: - output.hex_ids_.push_back(hex_id++); - break; - default: - SLIC_ERROR("unsupported element type"); - break; - } + switch (mesh.SpaceDimension()) { + case 2: + return Domain::ofElements(mesh, [](std::vector, int) { return true; }); + break; + case 3: + return Domain::ofElements(mesh, [](std::vector, int) { return true; }); + break; + default: + SLIC_ERROR("In valid spatial dimension. Domains may only be created on 2D or 3D meshes."); + exit(-1); } - - return output; } Domain EntireBoundary(const mfem::Mesh& mesh) { - Domain output{mesh, mesh.SpaceDimension() - 1, Domain::Type::BoundaryElements}; - - int edge_id = 0; - int tri_id = 0; - int quad_id = 0; - - for (int f = 0; f < mesh.GetNumFaces(); f++) { - // discard faces with the wrong type - if (mesh.GetFaceInformation(f).IsInterior()) continue; - - auto geom = mesh.GetFaceGeometry(f); - - switch (geom) { - case mfem::Geometry::SEGMENT: - output.edge_ids_.push_back(edge_id++); - break; - case mfem::Geometry::TRIANGLE: - output.tri_ids_.push_back(tri_id++); - break; - case mfem::Geometry::SQUARE: - output.quad_ids_.push_back(quad_id++); - break; - default: - SLIC_ERROR("unsupported element type"); - break; - } + switch (mesh.SpaceDimension()) { + case 2: + return Domain::ofBoundaryElements(mesh, [](std::vector, int) { return true; }); + break; + case 3: + return Domain::ofBoundaryElements(mesh, [](std::vector, int) { return true; }); + break; + default: + SLIC_ERROR("In valid spatial dimension. Domains may only be created on 2D or 3D meshes."); + exit(-1); } - - return output; } /// @cond @@ -553,22 +501,28 @@ Domain set_operation(set_op op, const Domain& a, const Domain& b) Domain output{a.mesh_, a.dim_}; - if (output.dim_ == 0) { - output.vertex_ids_ = set_operation(op, a.vertex_ids_, b.vertex_ids_); - } + using Ids = std::vector; + auto apply_set_op = [&op](const Ids& x, const Ids& y) { return set_operation(op, x, y); }; + + 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); + auto output_mfem_ids = apply_set_op(a_mfem_ids, b_mfem_ids); + output.addElements(output_ids, output_mfem_ids, g); + }; if (output.dim_ == 1) { - output.edge_ids_ = set_operation(op, a.edge_ids_, b.edge_ids_); + fill_output_lists(a.edge_ids_, a.mfem_edge_ids_, b.edge_ids_, b.mfem_edge_ids_, mfem::Geometry::SEGMENT); } if (output.dim_ == 2) { - output.tri_ids_ = set_operation(op, a.tri_ids_, b.tri_ids_); - output.quad_ids_ = set_operation(op, a.quad_ids_, b.quad_ids_); + fill_output_lists(a.tri_ids_, a.mfem_tri_ids_, b.tri_ids_, b.mfem_tri_ids_, mfem::Geometry::TRIANGLE); + fill_output_lists(a.quad_ids_, a.mfem_quad_ids_, b.quad_ids_, b.mfem_quad_ids_, mfem::Geometry::SQUARE); } if (output.dim_ == 3) { - output.tet_ids_ = set_operation(op, a.tet_ids_, b.tet_ids_); - output.hex_ids_ = set_operation(op, a.hex_ids_, b.hex_ids_); + fill_output_lists(a.tet_ids_, a.mfem_tet_ids_, b.tet_ids_, b.mfem_tet_ids_, mfem::Geometry::TETRAHEDRON); + fill_output_lists(a.hex_ids_, a.mfem_hex_ids_, b.hex_ids_, b.mfem_hex_ids_, mfem::Geometry::CUBE); } return output; diff --git a/src/serac/numerics/functional/domain.hpp b/src/serac/numerics/functional/domain.hpp index ad7c16911..7ecfdcfae 100644 --- a/src/serac/numerics/functional/domain.hpp +++ b/src/serac/numerics/functional/domain.hpp @@ -37,14 +37,38 @@ struct Domain { /// @brief whether the elements in this domain are on the boundary or not Type type_; - /// note: only lists with appropriate dimension (see dim_) will be populated - /// for example, a 2D Domain may have `tri_ids_` and `quad_ids_` non-nonempty, - /// but all other lists will be empty + ///@{ + /// @name ElementIds + /// Indices of elements contained in the domain. + /// 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 + /// hold indices into the "E-vector" of the appropriate geometry. These are + /// used primarily for identifying elements in the domain for participation + /// in integrals. /// - /// these lists hold indices into the "E-vector" of the appropriate geometry + /// The second set, (mfem_edge_ids_, mfem_tri_ids_, ...), gives the ids of + /// elements in this domain in the global mfem::Mesh data structure. These + /// maps are needed to find the dofs that live on a Domain. /// + /// Instances of Domain are meant to be homogeneous: only lists with + /// appropriate dimension (see dim_) will be populated by the factory + /// functions. For example, a 2D Domain may have `tri_ids_` and `quad_ids_` + /// non-empty, but all other lists will be empty. + /// + /// @note For every entry in the first group (say, edge_ids_), there should + /// be a corresponding entry into the second group (mfem_edge_ids_). This + /// is an intended invariant of the class, but it's not enforced by the data + /// structures. Prefer to use the factory methods (eg, \ref ofElements(...)) + /// to populate these lists automatically, as they repsect this invariant and + /// are tested. Otherwise, use the \ref addElements(...) or addElements(...) + /// 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 vertex_ids_; std::vector edge_ids_; std::vector tri_ids_; std::vector quad_ids_; @@ -58,6 +82,7 @@ struct Domain { 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) {} /** @@ -120,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_; @@ -132,6 +156,21 @@ struct Domain { /// @brief get mfem degree of freedom list for a given FiniteElementSpace mfem::Array dof_list(mfem::FiniteElementSpace* fes) const; + + /// @brief Add an element to the domain + /// + /// This is meant for internal use on the class. Prefer to use the factory + /// methods (ofElements, ofBoundaryElements, etc) to create domains and + /// thereby populate the element lists. + void addElement(int geom_id, int elem_id, mfem::Geometry::Type element_geometry); + + /// @brief Add a batch of elements to the domain + /// + /// This is meant for internal use on the class. Prefer to use the factory + /// methods (ofElements, ofBoundaryElements, etc) to create domains and + /// thereby populate the element lists. + void addElements(const std::vector& geom_id, const std::vector& elem_id, + mfem::Geometry::Type element_geometry); }; /// @brief constructs a domain from all the elements in a mesh diff --git a/src/serac/numerics/functional/tests/domain_tests.cpp b/src/serac/numerics/functional/tests/domain_tests.cpp index dfd1b0572..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) { { @@ -336,6 +276,134 @@ TEST(domain, of_elements) } } +TEST(domain, entireDomain2d) +{ + constexpr int dim = 2; + constexpr int p = 1; + auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); + + Domain d0 = EntireDomain(mesh); + + EXPECT_EQ(d0.dim_, 2); + EXPECT_EQ(d0.tri_ids_.size(), 2); + EXPECT_EQ(d0.quad_ids_.size(), 4); + + auto fec = mfem::H1_FECollection(p, dim); + auto fes = mfem::FiniteElementSpace(&mesh, &fec); + + mfem::Array dof_indices = d0.dof_list(&fes); + EXPECT_EQ(dof_indices.Size(), 8); +} + +TEST(domain, entireDomain3d) +{ + constexpr int dim = 3; + constexpr int p = 1; + auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); + + Domain d0 = EntireDomain(mesh); + + EXPECT_EQ(d0.dim_, 3); + EXPECT_EQ(d0.tet_ids_.size(), 12); + EXPECT_EQ(d0.hex_ids_.size(), 7); + + auto fec = mfem::H1_FECollection(p, dim); + auto fes = mfem::FiniteElementSpace(&mesh, &fec); + + mfem::Array dof_indices = d0.dof_list(&fes); + EXPECT_EQ(dof_indices.Size(), 25); +} + +TEST(domain, of2dElementsFindsDofs) +{ + constexpr int dim = 2; + constexpr int p = 2; + auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); + + auto fec = mfem::H1_FECollection(p, dim); + auto fes = mfem::FiniteElementSpace(&mesh, &fec); + + auto find_element_0 = [](std::vector vertices, int /* attr */) { + auto centroid = average(vertices); + return (centroid[0] < 0.5) && (centroid[1] < 0.25); + }; + + Domain d0 = Domain::ofElements(mesh, find_element_0); + + mfem::Array dof_indices = d0.dof_list(&fes); + + EXPECT_EQ(dof_indices.Size(), 9); + + /////////////////////////////////////// + + auto find_element_4 = [](std::vector vertices, int) { + auto centroid = average(vertices); + tensor target{{0.533, 0.424}}; + return norm(centroid - target) < 1e-2; + }; + Domain d1 = Domain::ofElements(mesh, find_element_4); + + Domain elements_0_and_4 = d0 | d1; + + dof_indices = elements_0_and_4.dof_list(&fes); + EXPECT_EQ(dof_indices.Size(), 12); + + /////////////////////////////////////// + + Domain d2 = EntireDomain(mesh) - elements_0_and_4; + + dof_indices = d2.dof_list(&fes); + + EXPECT_EQ(dof_indices.Size(), 22); +} + +TEST(domain, of3dElementsFindsDofs) +{ + constexpr int dim = 3; + constexpr int p = 2; + auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); + + auto fec = mfem::H1_FECollection(p, dim); + auto fes = mfem::FiniteElementSpace(&mesh, &fec); + + auto find_element_0 = [](std::vector vertices, int /* attr */) { + auto centroid = average(vertices); + vec3 target{{3.275, 0.7, 1.225}}; + return norm(centroid - target) < 1e-2; + }; + + Domain d0 = Domain::ofElements(mesh, find_element_0); + + mfem::Array dof_indices = d0.dof_list(&fes); + + // element 0 is a P2 tetrahedron, so it should have 10 dofs + EXPECT_EQ(dof_indices.Size(), 10); + + /////////////////////////////////////// + + auto find_element_1 = [](std::vector vertices, int) { + auto centroid = average(vertices); + vec3 target{{3.275, 1.2, 0.725}}; + return norm(centroid - target) < 1e-2; + }; + Domain d1 = Domain::ofElements(mesh, find_element_1); + + Domain elements_0_and_1 = d0 | d1; + + dof_indices = elements_0_and_1.dof_list(&fes); + + // Elements 0 and 1 are P2 tets that share one face -> 14 dofs + EXPECT_EQ(dof_indices.Size(), 14); + + ///////////////////////////////////////// + + Domain d2 = EntireDomain(mesh) - elements_0_and_1; + + dof_indices = d2.dof_list(&fes); + + EXPECT_EQ(dof_indices.Size(), 113); +} + int main(int argc, char* argv[]) { int num_procs, myid;