diff --git a/src/serac/numerics/functional/domain.hpp b/src/serac/numerics/functional/domain.hpp index 7ecfdcfae..d4c8bfb7f 100644 --- a/src/serac/numerics/functional/domain.hpp +++ b/src/serac/numerics/functional/domain.hpp @@ -18,7 +18,8 @@ namespace serac { * * This region can be an entire mesh or some subset of its elements */ -struct Domain { +class Domain { +public: /// @brief enum describing what kind of elements are included in a Domain enum Type { @@ -29,58 +30,13 @@ struct Domain { static constexpr int num_types = 2; ///< the number of entries in the Type enum /// @brief the underyling mesh for this domain - const mfem::Mesh& mesh_; + const mfem::Mesh& mesh() const { return mesh_; }; /// @brief the geometric dimension of the domain - int dim_; + int dim() const { return dim_; }; /// @brief whether the elements in this domain are on the boundary or not - Type type_; - - ///@{ - /// @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. - /// - /// 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 edge_ids_; - std::vector tri_ids_; - std::vector quad_ids_; - std::vector tet_ids_; - std::vector hex_ids_; - - std::vector mfem_edge_ids_; - std::vector mfem_tri_ids_; - std::vector mfem_quad_ids_; - std::vector mfem_tet_ids_; - std::vector mfem_hex_ids_; - /// @endcond + Type type() const { return type_; }; /// @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,6 +99,13 @@ struct Domain { static Domain ofBoundaryElements(const mfem::Mesh& mesh, std::function, int)> func); /// @brief get elements by geometry type + /// + /// Returned list holds the index of each element in this Domain in the set + /// of all elements of like geometry in the mesh. + /// For example, if get(mfem::Geometry::SEGMENT)[0] is 5, then element 0 in + /// this domain is element 5 in the grouping of all edges in the mesh. In + /// other words, the returned list holds indices into the "E-vector" of the + /// appropriate geometry. const std::vector& get(mfem::Geometry::Type geom) const { if (geom == mfem::Geometry::SEGMENT) return edge_ids_; @@ -157,20 +120,87 @@ 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. + /** @brief Add an element to the domain + * + * Prefer to use the factory methods (ofElements, ofBoundaryElements, etc) + * to create domains and thereby populate the element lists. Does not + * check validity of inputs. + */ 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. + /** @brief Add a batch of elements to the domain + * + * Prefer to use the factory methods (ofElements, ofBoundaryElements, etc) + * to create domains and thereby populate the element lists. Does not + * check validity of inputs. + */ void addElements(const std::vector& geom_id, const std::vector& elem_id, mfem::Geometry::Type element_geometry); + +private: + /// @brief the underyling mesh for this domain + const mfem::Mesh& mesh_; + + /// @brief the geometric dimension of the domain + int dim_; + + /// @brief whether the elements in this domain are on the boundary or not + Type type_; + + ///@{ + /// @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. + /// + /// 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 edge_ids_; + std::vector tri_ids_; + std::vector quad_ids_; + std::vector tet_ids_; + std::vector hex_ids_; + + std::vector mfem_edge_ids_; + std::vector mfem_tri_ids_; + std::vector mfem_quad_ids_; + std::vector mfem_tet_ids_; + std::vector mfem_hex_ids_; + /// @endcond + + using c_iter = std::vector::const_iterator; + using b_iter = std::back_insert_iterator>; + using set_op = std::function; + + friend Domain set_operation(set_op op, const Domain& a, const Domain& b); + + friend Domain operator|(const Domain& a, const Domain& b); + friend Domain operator&(const Domain& a, const Domain& b); + friend Domain operator-(const Domain& a, const Domain& b); }; /// @brief constructs a domain from all the elements in a mesh @@ -192,7 +222,7 @@ Domain operator-(const Domain& a, const Domain& b); template inline auto by_attr(int value) { - return [value](std::vector >, int attr) { return attr == value; }; + return [value](std::vector>, int attr) { return attr == value; }; } } // namespace serac diff --git a/src/serac/numerics/functional/functional.hpp b/src/serac/numerics/functional/functional.hpp index 79311b68c..45dcdac78 100644 --- a/src/serac/numerics/functional/functional.hpp +++ b/src/serac/numerics/functional/functional.hpp @@ -304,12 +304,12 @@ class Functional { void AddDomainIntegral(Dimension, DependsOn, const lambda& integrand, const Domain& domain, std::shared_ptr> qdata = NoQData) { - if (domain.mesh_.GetNE() == 0) return; + if (domain.mesh().GetNE() == 0) return; - SLIC_ERROR_ROOT_IF(dim != domain.mesh_.Dimension(), "invalid mesh dimension for domain integral"); + SLIC_ERROR_ROOT_IF(dim != domain.mesh().Dimension(), "invalid mesh dimension for domain integral"); - check_for_unsupported_elements(domain.mesh_); - check_for_missing_nodal_gridfunc(domain.mesh_); + check_for_unsupported_elements(domain.mesh()); + check_for_missing_nodal_gridfunc(domain.mesh()); using signature = test(decltype(serac::type(trial_spaces))...); integrals_.push_back( @@ -342,12 +342,12 @@ class Functional { template void AddBoundaryIntegral(Dimension, DependsOn, const lambda& integrand, const Domain& domain) { - auto num_bdr_elements = domain.mesh_.GetNBE(); + auto num_bdr_elements = domain.mesh().GetNBE(); if (num_bdr_elements == 0) return; - SLIC_ERROR_ROOT_IF(dim != domain.dim_, "invalid domain of integration for boundary integral"); + SLIC_ERROR_ROOT_IF(dim != domain.dim(), "invalid domain of integration for boundary integral"); - check_for_missing_nodal_gridfunc(domain.mesh_); + check_for_missing_nodal_gridfunc(domain.mesh()); using signature = test(decltype(serac::type(trial_spaces))...); integrals_.push_back(MakeBoundaryIntegral(domain, integrand, std::vector{args...})); @@ -414,7 +414,7 @@ class Functional { bool already_computed[Domain::num_types]{}; // default initializes to `false` for (auto& integral : integrals_) { - auto type = integral.domain_.type_; + auto type = integral.domain_.type(); if (!already_computed[type]) { G_trial_[type][which].Gather(input_L_[which], input_E_[type][which]); @@ -460,7 +460,7 @@ class Functional { bool already_computed[Domain::num_types][num_trial_spaces]{}; // default initializes to `false` for (auto& integral : integrals_) { - auto type = integral.domain_.type_; + auto type = integral.domain_.type(); for (auto i : integral.active_trial_spaces_) { if (!already_computed[type][i]) { @@ -588,9 +588,9 @@ class Functional { std::map> element_gradients[Domain::num_types]; for (auto& integral : form_.integrals_) { - auto& K_elem = element_gradients[integral.domain_.type_]; - auto& test_restrictions = form_.G_test_[integral.domain_.type_].restrictions; - auto& trial_restrictions = form_.G_trial_[integral.domain_.type_][which_argument].restrictions; + auto& K_elem = element_gradients[integral.domain_.type()]; + auto& test_restrictions = form_.G_test_[integral.domain_.type()].restrictions; + auto& trial_restrictions = form_.G_trial_[integral.domain_.type()][which_argument].restrictions; if (K_elem.empty()) { for (auto& [geom, test_restriction] : test_restrictions) { diff --git a/src/serac/numerics/functional/functional_qoi.inl b/src/serac/numerics/functional/functional_qoi.inl index 2c0cea577..c7dfec37d 100644 --- a/src/serac/numerics/functional/functional_qoi.inl +++ b/src/serac/numerics/functional/functional_qoi.inl @@ -177,12 +177,12 @@ public: void AddDomainIntegral(Dimension, DependsOn, const lambda& integrand, Domain& domain, std::shared_ptr> qdata = NoQData) { - if (domain.mesh_.GetNE() == 0) return; + if (domain.mesh().GetNE() == 0) return; - SLIC_ERROR_ROOT_IF(dim != domain.mesh_.Dimension(), "invalid mesh dimension for domain integral"); + SLIC_ERROR_ROOT_IF(dim != domain.mesh().Dimension(), "invalid mesh dimension for domain integral"); - check_for_unsupported_elements(domain.mesh_); - check_for_missing_nodal_gridfunc(domain.mesh_); + check_for_unsupported_elements(domain.mesh()); + check_for_missing_nodal_gridfunc(domain.mesh()); using signature = test(decltype(serac::type(trial_spaces))...); integrals_.push_back( @@ -217,12 +217,12 @@ public: template void AddBoundaryIntegral(Dimension, DependsOn, const lambda& integrand, const Domain& domain) { - auto num_bdr_elements = domain.mesh_.GetNBE(); + auto num_bdr_elements = domain.mesh().GetNBE(); if (num_bdr_elements == 0) return; - SLIC_ERROR_ROOT_IF(dim != domain.dim_, "invalid domain of integration for boundary integral"); + SLIC_ERROR_ROOT_IF(dim != domain.dim(), "invalid domain of integration for boundary integral"); - check_for_missing_nodal_gridfunc(domain.mesh_); + check_for_missing_nodal_gridfunc(domain.mesh()); using signature = test(decltype(serac::type(trial_spaces))...); integrals_.push_back(MakeBoundaryIntegral(domain, integrand, std::vector{args...})); @@ -290,7 +290,7 @@ public: bool already_computed[Domain::num_types]{}; // default initializes to `false` for (auto& integral : integrals_) { - auto type = integral.domain_.type_; + auto type = integral.domain_.type(); if (!already_computed[type]) { G_trial_[type][which].Gather(input_L_[which], input_E_[type][which]); @@ -336,7 +336,7 @@ public: bool already_computed[Domain::num_types][num_trial_spaces]{}; // default initializes to `false` for (auto& integral : integrals_) { - auto type = integral.domain_.type_; + auto type = integral.domain_.type(); for (auto i : integral.active_trial_spaces_) { if (!already_computed[type][i]) { @@ -430,8 +430,8 @@ private: std::map> element_gradients[Domain::num_types]; for (auto& integral : form_.integrals_) { - auto& K_elem = element_gradients[integral.domain_.type_]; - auto& trial_restrictions = form_.G_trial_[integral.domain_.type_][which_argument].restrictions; + auto& K_elem = element_gradients[integral.domain_.type()]; + auto& trial_restrictions = form_.G_trial_[integral.domain_.type()][which_argument].restrictions; if (K_elem.empty()) { for (auto& [geom, trial_restriction] : trial_restrictions) { diff --git a/src/serac/numerics/functional/geometric_factors.cpp b/src/serac/numerics/functional/geometric_factors.cpp index 7e8897fa5..4ff7147d3 100644 --- a/src/serac/numerics/functional/geometric_factors.cpp +++ b/src/serac/numerics/functional/geometric_factors.cpp @@ -61,7 +61,7 @@ void compute_geometric_factors(mfem::Vector& positions_q, mfem::Vector& jacobian GeometricFactors::GeometricFactors(const Domain& d, int q, mfem::Geometry::Type g) { - auto* nodes = d.mesh_.GetNodes(); + auto* nodes = d.mesh().GetNodes(); auto* fes = nodes->FESpace(); auto restriction = serac::ElementRestriction(fes, g); @@ -71,14 +71,11 @@ GeometricFactors::GeometricFactors(const Domain& d, int q, mfem::Geometry::Type // assumes all elements are the same order int p = fes->GetElementOrder(0); - int spatial_dim = d.mesh_.SpaceDimension(); + int spatial_dim = d.mesh().SpaceDimension(); int geometry_dim = dimension_of(g); int qpts_per_elem = num_quadrature_points(g, q); - if (g == mfem::Geometry::TRIANGLE) elements = d.tri_ids_; - if (g == mfem::Geometry::SQUARE) elements = d.quad_ids_; - if (g == mfem::Geometry::TETRAHEDRON) elements = d.tet_ids_; - if (g == mfem::Geometry::CUBE) elements = d.hex_ids_; + elements = d.get(g); num_elements = elements.size(); @@ -139,7 +136,7 @@ GeometricFactors::GeometricFactors(const Domain& d, int q, mfem::Geometry::Type GeometricFactors::GeometricFactors(const Domain& d, int q, mfem::Geometry::Type g, FaceType type) { - auto* nodes = d.mesh_.GetNodes(); + auto* nodes = d.mesh().GetNodes(); auto* fes = nodes->FESpace(); auto restriction = serac::ElementRestriction(fes, g, type); @@ -149,7 +146,7 @@ GeometricFactors::GeometricFactors(const Domain& d, int q, mfem::Geometry::Type // assumes all elements are the same order int p = fes->GetElementOrder(0); - int spatial_dim = d.mesh_.SpaceDimension(); + int spatial_dim = d.mesh().SpaceDimension(); int geometry_dim = dimension_of(g); int qpts_per_elem = num_quadrature_points(g, q); diff --git a/src/serac/numerics/functional/integral.hpp b/src/serac/numerics/functional/integral.hpp index 337118a7d..df10c293d 100644 --- a/src/serac/numerics/functional/integral.hpp +++ b/src/serac/numerics/functional/integral.hpp @@ -242,7 +242,7 @@ Integral MakeDomainIntegral(const Domain& domain, const lambda_type& qf, { FunctionSignature signature; - SLIC_ERROR_IF(domain.type_ != Domain::Type::Elements, "Error: trying to evaluate a domain integral over a boundary"); + SLIC_ERROR_IF(domain.type() != Domain::Type::Elements, "Error: trying to evaluate a domain integral over a boundary"); Integral integral(domain, argument_indices); @@ -329,7 +329,7 @@ Integral MakeBoundaryIntegral(const Domain& domain, const lambda_type& qf, std:: { FunctionSignature signature; - SLIC_ERROR_IF(domain.type_ != Domain::Type::BoundaryElements, + SLIC_ERROR_IF(domain.type() != Domain::Type::BoundaryElements, "Error: trying to evaluate a boundary integral over a non-boundary domain of integration"); Integral integral(domain, argument_indices); diff --git a/src/serac/numerics/functional/tests/domain_tests.cpp b/src/serac/numerics/functional/tests/domain_tests.cpp index 7893107a9..538960e9c 100644 --- a/src/serac/numerics/functional/tests/domain_tests.cpp +++ b/src/serac/numerics/functional/tests/domain_tests.cpp @@ -34,26 +34,30 @@ tensor average(std::vector >& positions) TEST(domain, of_edges) { { - auto mesh = import_mesh("onehex.mesh"); - Domain d0 = Domain::ofEdges(mesh, std::function([](std::vector x) { + auto mesh = import_mesh("onehex.mesh"); + Domain d0 = Domain::ofEdges(mesh, std::function([](std::vector x) { return (0.5 * (x[0][0] + x[1][0])) < 0.25; // x coordinate of edge midpoint })); - EXPECT_EQ(d0.edge_ids_.size(), 4); - EXPECT_EQ(d0.dim_, 1); + auto d0_edges = d0.get(mfem::Geometry::SEGMENT); + EXPECT_EQ(d0_edges.size(), 4); + EXPECT_EQ(d0.dim(), 1); - Domain d1 = Domain::ofEdges(mesh, std::function([](std::vector x) { + Domain d1 = Domain::ofEdges(mesh, std::function([](std::vector x) { return (0.5 * (x[0][1] + x[1][1])) < 0.25; // y coordinate of edge midpoint })); - EXPECT_EQ(d1.edge_ids_.size(), 4); - EXPECT_EQ(d1.dim_, 1); + auto d1_edges = d1.get(mfem::Geometry::SEGMENT); + EXPECT_EQ(d1_edges.size(), 4); + EXPECT_EQ(d1.dim(), 1); - Domain d2 = d0 | d1; - EXPECT_EQ(d2.edge_ids_.size(), 7); - EXPECT_EQ(d2.dim_, 1); + Domain d2 = d0 | d1; + auto d2_edges = d2.get(mfem::Geometry::SEGMENT); + EXPECT_EQ(d2_edges.size(), 7); + EXPECT_EQ(d2.dim(), 1); - Domain d3 = d0 & d1; - EXPECT_EQ(d3.edge_ids_.size(), 1); - EXPECT_EQ(d3.dim_, 1); + Domain d3 = d0 & d1; + auto d3_edges = d3.get(mfem::Geometry::SEGMENT); + EXPECT_EQ(d3_edges.size(), 1); + EXPECT_EQ(d3.dim(), 1); // note: by_attr doesn't apply to edge sets in 3D, since // mfem doesn't have the notion of edge attributes @@ -61,26 +65,30 @@ TEST(domain, of_edges) } { - auto mesh = import_mesh("onetet.mesh"); - Domain d0 = Domain::ofEdges(mesh, std::function([](std::vector x) { + auto mesh = import_mesh("onetet.mesh"); + Domain d0 = Domain::ofEdges(mesh, std::function([](std::vector x) { return (0.5 * (x[0][0] + x[1][0])) < 0.25; // x coordinate of edge midpoint })); - EXPECT_EQ(d0.edge_ids_.size(), 3); - EXPECT_EQ(d0.dim_, 1); + auto d0_edges = d0.get(mfem::Geometry::SEGMENT); + EXPECT_EQ(d0_edges.size(), 3); + EXPECT_EQ(d0.dim(), 1); - Domain d1 = Domain::ofEdges(mesh, std::function([](std::vector x) { + Domain d1 = Domain::ofEdges(mesh, std::function([](std::vector x) { return (0.5 * (x[0][1] + x[1][1])) < 0.25; // y coordinate of edge midpoint })); - EXPECT_EQ(d1.edge_ids_.size(), 3); - EXPECT_EQ(d1.dim_, 1); + auto d1_edges = d1.get(mfem::Geometry::SEGMENT); + EXPECT_EQ(d1_edges.size(), 3); + EXPECT_EQ(d1.dim(), 1); - Domain d2 = d0 | d1; - EXPECT_EQ(d2.edge_ids_.size(), 5); - EXPECT_EQ(d2.dim_, 1); + Domain d2 = d0 | d1; + auto d2_edges = d2.get(mfem::Geometry::SEGMENT); + EXPECT_EQ(d2_edges.size(), 5); + EXPECT_EQ(d2.dim(), 1); - Domain d3 = d0 & d1; - EXPECT_EQ(d3.edge_ids_.size(), 1); - EXPECT_EQ(d3.dim_, 1); + Domain d3 = d0 & d1; + auto d3_edges = d3.get(mfem::Geometry::SEGMENT); + EXPECT_EQ(d3_edges.size(), 1); + EXPECT_EQ(d3.dim(), 1); // note: by_attr doesn't apply to edge sets in 3D, since // mfem doesn't have the notion of edge attributes @@ -91,121 +99,140 @@ TEST(domain, of_edges) constexpr int dim = 2; auto mesh = import_mesh("beam-quad.mesh"); mesh.FinalizeQuadMesh(true); - Domain d0 = Domain::ofEdges(mesh, std::function([](std::vector x, int /* bdr_attr */) { + Domain d0 = Domain::ofEdges(mesh, std::function([](std::vector x, int /* bdr_attr */) { return (0.5 * (x[0][0] + x[1][0])) < 0.25; // x coordinate of edge midpoint })); - EXPECT_EQ(d0.edge_ids_.size(), 1); - EXPECT_EQ(d0.dim_, 1); + auto d0_edges = d0.get(mfem::Geometry::SEGMENT); + EXPECT_EQ(d0_edges.size(), 1); + EXPECT_EQ(d0.dim(), 1); - Domain d1 = Domain::ofEdges(mesh, std::function([](std::vector x, int /* bdr_attr */) { + Domain d1 = Domain::ofEdges(mesh, std::function([](std::vector x, int /* bdr_attr */) { return (0.5 * (x[0][1] + x[1][1])) < 0.25; // y coordinate of edge midpoint })); - EXPECT_EQ(d1.edge_ids_.size(), 8); - EXPECT_EQ(d1.dim_, 1); + auto d1_edges = d1.get(mfem::Geometry::SEGMENT); + EXPECT_EQ(d1_edges.size(), 8); + EXPECT_EQ(d1.dim(), 1); - Domain d2 = d0 | d1; - EXPECT_EQ(d2.edge_ids_.size(), 9); - EXPECT_EQ(d2.dim_, 1); + Domain d2 = d0 | d1; + auto d2_edges = d2.get(mfem::Geometry::SEGMENT); + EXPECT_EQ(d2_edges.size(), 9); + EXPECT_EQ(d2.dim(), 1); - Domain d3 = d0 & d1; - EXPECT_EQ(d3.edge_ids_.size(), 0); - EXPECT_EQ(d3.dim_, 1); + Domain d3 = d0 & d1; + auto d3_edges = d3.get(mfem::Geometry::SEGMENT); + EXPECT_EQ(d3_edges.size(), 0); + EXPECT_EQ(d3.dim(), 1); // check that by_attr compiles Domain d4 = Domain::ofEdges(mesh, by_attr(3)); - Domain d5 = Domain::ofBoundaryElements(mesh, [](std::vector, int) { return true; }); - EXPECT_EQ(d5.edge_ids_.size(), 18); // 1x8 row of quads has 18 boundary edges + Domain d5 = Domain::ofBoundaryElements(mesh, [](std::vector, int) { return true; }); + auto d5_edges = d5.get(mfem::Geometry::SEGMENT); + EXPECT_EQ(d5_edges.size(), 18); // 1x8 row of quads has 18 boundary edges } } TEST(domain, of_faces) { { - constexpr int dim = 3; - auto mesh = import_mesh("onehex.mesh"); - Domain d0 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { + constexpr int dim = 3; + auto mesh = import_mesh("onehex.mesh"); + Domain d0 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { return average(vertices)[0] < 0.25; // x coordinate of face center })); - EXPECT_EQ(d0.quad_ids_.size(), 1); - EXPECT_EQ(d0.dim_, 2); + auto d0_quads = d0.get(mfem::Geometry::SQUARE); + EXPECT_EQ(d0_quads.size(), 1); + EXPECT_EQ(d0.dim(), 2); - Domain d1 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { + Domain d1 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { return average(vertices)[1] < 0.25; // y coordinate of face center })); - EXPECT_EQ(d1.quad_ids_.size(), 1); - EXPECT_EQ(d1.dim_, 2); + auto d1_quads = d1.get(mfem::Geometry::SQUARE); + EXPECT_EQ(d1_quads.size(), 1); + EXPECT_EQ(d1.dim(), 2); - Domain d2 = d0 | d1; - EXPECT_EQ(d2.quad_ids_.size(), 2); - EXPECT_EQ(d2.dim_, 2); + Domain d2 = d0 | d1; + auto d2_quads = d2.get(mfem::Geometry::SQUARE); + EXPECT_EQ(d2_quads.size(), 2); + EXPECT_EQ(d2.dim(), 2); - Domain d3 = d0 & d1; - EXPECT_EQ(d3.quad_ids_.size(), 0); - EXPECT_EQ(d3.dim_, 2); + Domain d3 = d0 & d1; + auto d3_quads = d3.get(mfem::Geometry::SQUARE); + EXPECT_EQ(d3_quads.size(), 0); + EXPECT_EQ(d3.dim(), 2); // check that by_attr compiles Domain d4 = Domain::ofFaces(mesh, by_attr(3)); - Domain d5 = Domain::ofBoundaryElements(mesh, [](std::vector, int) { return true; }); - EXPECT_EQ(d5.quad_ids_.size(), 6); + Domain d5 = Domain::ofBoundaryElements(mesh, [](std::vector, int) { return true; }); + auto d5_quads = d5.get(mfem::Geometry::SQUARE); + EXPECT_EQ(d5_quads.size(), 6); } { - constexpr int dim = 3; - auto mesh = import_mesh("onetet.mesh"); - Domain d0 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /* bdr_attr */) { + constexpr int dim = 3; + auto mesh = import_mesh("onetet.mesh"); + Domain d0 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /* bdr_attr */) { // accept face if it contains a vertex whose x coordinate is less than 0.1 for (auto v : vertices) { if (v[0] < 0.1) return true; } return false; })); - EXPECT_EQ(d0.tri_ids_.size(), 4); - EXPECT_EQ(d0.dim_, 2); + auto d0_tris = d0.get(mfem::Geometry::TRIANGLE); + EXPECT_EQ(d0_tris.size(), 4); + EXPECT_EQ(d0.dim(), 2); Domain d1 = Domain::ofFaces( mesh, std::function([](std::vector x, int /* bdr_attr */) { return average(x)[1] < 0.1; })); - EXPECT_EQ(d1.tri_ids_.size(), 1); - EXPECT_EQ(d1.dim_, 2); + auto d1_tris = d1.get(mfem::Geometry::TRIANGLE); + EXPECT_EQ(d1_tris.size(), 1); + EXPECT_EQ(d1.dim(), 2); - Domain d2 = d0 | d1; - EXPECT_EQ(d2.tri_ids_.size(), 4); - EXPECT_EQ(d2.dim_, 2); + Domain d2 = d0 | d1; + auto d2_tris = d2.get(mfem::Geometry::TRIANGLE); + EXPECT_EQ(d2_tris.size(), 4); + EXPECT_EQ(d2.dim(), 2); - Domain d3 = d0 & d1; - EXPECT_EQ(d3.tri_ids_.size(), 1); - EXPECT_EQ(d3.dim_, 2); + Domain d3 = d0 & d1; + auto d3_tris = d3.get(mfem::Geometry::TRIANGLE); + EXPECT_EQ(d3_tris.size(), 1); + EXPECT_EQ(d3.dim(), 2); // check that by_attr compiles Domain d4 = Domain::ofFaces(mesh, by_attr(3)); - Domain d5 = Domain::ofBoundaryElements(mesh, [](std::vector, int) { return true; }); - EXPECT_EQ(d5.tri_ids_.size(), 4); + Domain d5 = Domain::ofBoundaryElements(mesh, [](std::vector, int) { return true; }); + auto d5_tris = d5.get(mfem::Geometry::TRIANGLE); + EXPECT_EQ(d5_tris.size(), 4); } { - constexpr int dim = 2; - auto mesh = import_mesh("beam-quad.mesh"); - Domain d0 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /* attr */) { + constexpr int dim = 2; + auto mesh = import_mesh("beam-quad.mesh"); + Domain d0 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[0] < 2.25; // x coordinate of face center })); - EXPECT_EQ(d0.quad_ids_.size(), 2); - EXPECT_EQ(d0.dim_, 2); + auto d0_quads = d0.get(mfem::Geometry::SQUARE); + EXPECT_EQ(d0_quads.size(), 2); + EXPECT_EQ(d0.dim(), 2); - Domain d1 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /* attr */) { + Domain d1 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[1] < 0.55; // y coordinate of face center })); - EXPECT_EQ(d1.quad_ids_.size(), 8); - EXPECT_EQ(d1.dim_, 2); + auto d1_quads = d1.get(mfem::Geometry::SQUARE); + EXPECT_EQ(d1_quads.size(), 8); + EXPECT_EQ(d1.dim(), 2); - Domain d2 = d0 | d1; - EXPECT_EQ(d2.quad_ids_.size(), 8); - EXPECT_EQ(d2.dim_, 2); + Domain d2 = d0 | d1; + auto d2_quads = d2.get(mfem::Geometry::SQUARE); + EXPECT_EQ(d2_quads.size(), 8); + EXPECT_EQ(d2.dim(), 2); - Domain d3 = d0 & d1; - EXPECT_EQ(d3.quad_ids_.size(), 2); - EXPECT_EQ(d3.dim_, 2); + Domain d3 = d0 & d1; + auto d3_quads = d3.get(mfem::Geometry::SQUARE); + EXPECT_EQ(d3_quads.size(), 2); + EXPECT_EQ(d3.dim(), 2); // check that by_attr compiles Domain d4 = Domain::ofFaces(mesh, by_attr(3)); @@ -221,26 +248,26 @@ TEST(domain, of_elements) return average(vertices)[0] < 0.7; // x coordinate of face center })); - EXPECT_EQ(d0.tet_ids_.size(), 0); - EXPECT_EQ(d0.hex_ids_.size(), 1); - EXPECT_EQ(d0.dim_, 3); + EXPECT_EQ(d0.get(mfem::Geometry::TETRAHEDRON).size(), 0); + EXPECT_EQ(d0.get(mfem::Geometry::CUBE).size(), 1); + EXPECT_EQ(d0.dim(), 3); Domain d1 = Domain::ofElements(mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { return average(vertices)[1] < 0.75; // y coordinate of face center })); - EXPECT_EQ(d1.tet_ids_.size(), 6); - EXPECT_EQ(d1.hex_ids_.size(), 1); - EXPECT_EQ(d1.dim_, 3); + EXPECT_EQ(d1.get(mfem::Geometry::TETRAHEDRON).size(), 6); + EXPECT_EQ(d1.get(mfem::Geometry::CUBE).size(), 1); + EXPECT_EQ(d1.dim(), 3); Domain d2 = d0 | d1; - EXPECT_EQ(d2.tet_ids_.size(), 6); - EXPECT_EQ(d2.hex_ids_.size(), 2); - EXPECT_EQ(d2.dim_, 3); + EXPECT_EQ(d2.get(mfem::Geometry::TETRAHEDRON).size(), 6); + EXPECT_EQ(d2.get(mfem::Geometry::CUBE).size(), 2); + EXPECT_EQ(d2.dim(), 3); Domain d3 = d0 & d1; - EXPECT_EQ(d3.tet_ids_.size(), 0); - EXPECT_EQ(d3.hex_ids_.size(), 0); - EXPECT_EQ(d3.dim_, 3); + EXPECT_EQ(d3.get(mfem::Geometry::TETRAHEDRON).size(), 0); + EXPECT_EQ(d3.get(mfem::Geometry::CUBE).size(), 0); + EXPECT_EQ(d3.dim(), 3); // check that by_attr works Domain d4 = Domain::ofElements(mesh, by_attr(3)); @@ -251,25 +278,25 @@ TEST(domain, of_elements) auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); Domain d0 = Domain::ofElements( mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[0] < 0.45; })); - EXPECT_EQ(d0.tri_ids_.size(), 1); - EXPECT_EQ(d0.quad_ids_.size(), 1); - EXPECT_EQ(d0.dim_, 2); + EXPECT_EQ(d0.get(mfem::Geometry::TRIANGLE).size(), 1); + EXPECT_EQ(d0.get(mfem::Geometry::SQUARE).size(), 1); + EXPECT_EQ(d0.dim(), 2); Domain d1 = Domain::ofElements( mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[1] < 0.45; })); - EXPECT_EQ(d1.tri_ids_.size(), 1); - EXPECT_EQ(d1.quad_ids_.size(), 1); - EXPECT_EQ(d1.dim_, 2); + EXPECT_EQ(d1.get(mfem::Geometry::TRIANGLE).size(), 1); + EXPECT_EQ(d1.get(mfem::Geometry::SQUARE).size(), 1); + EXPECT_EQ(d1.dim(), 2); Domain d2 = d0 | d1; - EXPECT_EQ(d2.tri_ids_.size(), 2); - EXPECT_EQ(d2.quad_ids_.size(), 2); - EXPECT_EQ(d2.dim_, 2); + EXPECT_EQ(d2.get(mfem::Geometry::TRIANGLE).size(), 2); + EXPECT_EQ(d2.get(mfem::Geometry::SQUARE).size(), 2); + EXPECT_EQ(d2.dim(), 2); Domain d3 = d0 & d1; - EXPECT_EQ(d3.tri_ids_.size(), 0); - EXPECT_EQ(d3.quad_ids_.size(), 0); - EXPECT_EQ(d3.dim_, 2); + EXPECT_EQ(d3.get(mfem::Geometry::TRIANGLE).size(), 0); + EXPECT_EQ(d3.get(mfem::Geometry::SQUARE).size(), 0); + EXPECT_EQ(d3.dim(), 2); // check that by_attr compiles Domain d4 = Domain::ofElements(mesh, by_attr(3)); @@ -284,9 +311,9 @@ TEST(domain, entireDomain2d) Domain d0 = EntireDomain(mesh); - EXPECT_EQ(d0.dim_, 2); - EXPECT_EQ(d0.tri_ids_.size(), 2); - EXPECT_EQ(d0.quad_ids_.size(), 4); + EXPECT_EQ(d0.dim(), 2); + EXPECT_EQ(d0.get(mfem::Geometry::TRIANGLE).size(), 2); + EXPECT_EQ(d0.get(mfem::Geometry::SQUARE).size(), 4); auto fec = mfem::H1_FECollection(p, dim); auto fes = mfem::FiniteElementSpace(&mesh, &fec); @@ -303,9 +330,9 @@ TEST(domain, entireDomain3d) Domain d0 = EntireDomain(mesh); - EXPECT_EQ(d0.dim_, 3); - EXPECT_EQ(d0.tet_ids_.size(), 12); - EXPECT_EQ(d0.hex_ids_.size(), 7); + EXPECT_EQ(d0.dim(), 3); + EXPECT_EQ(d0.get(mfem::Geometry::TETRAHEDRON).size(), 12); + EXPECT_EQ(d0.get(mfem::Geometry::CUBE).size(), 7); auto fec = mfem::H1_FECollection(p, dim); auto fes = mfem::FiniteElementSpace(&mesh, &fec);