Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Hide data in Domain #1277

Open
wants to merge 6 commits into
base: develop
Choose a base branch
from
Open

Hide data in Domain #1277

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
d307e4e
Make element id vectors private
btalamini Nov 29, 2024
8efeb17
Make remaining member variables private
btalamini Nov 29, 2024
9c4085a
Refactor to eliminate unneeded class I introduced
btalamini Nov 29, 2024
8e43e94
Update docstrings
btalamini Nov 29, 2024
cae16eb
Make style
btalamini Nov 29, 2024
c38efa0
Fix doxygen error
btalamini Nov 29, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
150 changes: 90 additions & 60 deletions src/serac/numerics/functional/domain.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
{
Expand All @@ -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<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_edge_ids_;
std::vector<int> mfem_tri_ids_;
std::vector<int> mfem_quad_ids_;
std::vector<int> mfem_tet_ids_;
std::vector<int> 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) {}
Expand Down Expand Up @@ -143,6 +99,13 @@ struct Domain {
static Domain ofBoundaryElements(const mfem::Mesh& mesh, std::function<bool(std::vector<vec3>, 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<int>& get(mfem::Geometry::Type geom) const
{
if (geom == mfem::Geometry::SEGMENT) return edge_ids_;
Expand All @@ -157,20 +120,87 @@ struct Domain {
/// @brief get mfem degree of freedom list for a given FiniteElementSpace
mfem::Array<int> 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<int>& geom_id, const std::vector<int>& 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<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_edge_ids_;
std::vector<int> mfem_tri_ids_;
std::vector<int> mfem_quad_ids_;
std::vector<int> mfem_tet_ids_;
std::vector<int> mfem_hex_ids_;
/// @endcond

using c_iter = std::vector<int>::const_iterator;
using b_iter = std::back_insert_iterator<std::vector<int>>;
using set_op = std::function<b_iter(c_iter, c_iter, c_iter, c_iter, b_iter)>;

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
Expand All @@ -192,7 +222,7 @@ Domain operator-(const Domain& a, const Domain& b);
template <int dim>
inline auto by_attr(int value)
{
return [value](std::vector<tensor<double, dim> >, int attr) { return attr == value; };
return [value](std::vector<tensor<double, dim>>, int attr) { return attr == value; };
}

} // namespace serac
24 changes: 12 additions & 12 deletions src/serac/numerics/functional/functional.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -304,12 +304,12 @@ class Functional<test(trials...), exec> {
void AddDomainIntegral(Dimension<dim>, DependsOn<args...>, const lambda& integrand, const Domain& domain,
std::shared_ptr<QuadratureData<qpt_data_type>> 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<args>(trial_spaces))...);
integrals_.push_back(
Expand Down Expand Up @@ -342,12 +342,12 @@ class Functional<test(trials...), exec> {
template <int dim, int... args, typename lambda>
void AddBoundaryIntegral(Dimension<dim>, DependsOn<args...>, 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<args>(trial_spaces))...);
integrals_.push_back(MakeBoundaryIntegral<signature, Q, dim>(domain, integrand, std::vector<uint32_t>{args...}));
Expand Down Expand Up @@ -414,7 +414,7 @@ class Functional<test(trials...), exec> {
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]);
Expand Down Expand Up @@ -460,7 +460,7 @@ class Functional<test(trials...), exec> {
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]) {
Expand Down Expand Up @@ -588,9 +588,9 @@ class Functional<test(trials...), exec> {
std::map<mfem::Geometry::Type, ExecArray<double, 3, exec>> 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) {
Expand Down
22 changes: 11 additions & 11 deletions src/serac/numerics/functional/functional_qoi.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -177,12 +177,12 @@ public:
void AddDomainIntegral(Dimension<dim>, DependsOn<args...>, const lambda& integrand, Domain& domain,
std::shared_ptr<QuadratureData<qpt_data_type>> 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<args>(trial_spaces))...);
integrals_.push_back(
Expand Down Expand Up @@ -217,12 +217,12 @@ public:
template <int dim, int... args, typename lambda>
void AddBoundaryIntegral(Dimension<dim>, DependsOn<args...>, 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<args>(trial_spaces))...);
integrals_.push_back(MakeBoundaryIntegral<signature, Q, dim>(domain, integrand, std::vector<uint32_t>{args...}));
Expand Down Expand Up @@ -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]);
Expand Down Expand Up @@ -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]) {
Expand Down Expand Up @@ -430,8 +430,8 @@ private:
std::map<mfem::Geometry::Type, ExecArray<double, 3, exec>> 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) {
Expand Down
13 changes: 5 additions & 8 deletions src/serac/numerics/functional/geometric_factors.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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();

Expand Down Expand Up @@ -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);
Expand All @@ -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);

Expand Down
4 changes: 2 additions & 2 deletions src/serac/numerics/functional/integral.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -242,7 +242,7 @@ Integral MakeDomainIntegral(const Domain& domain, const lambda_type& qf,
{
FunctionSignature<s> 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);

Expand Down Expand Up @@ -329,7 +329,7 @@ Integral MakeBoundaryIntegral(const Domain& domain, const lambda_type& qf, std::
{
FunctionSignature<s> 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);
Expand Down
Loading