Skip to content

Fix acquisition of dof lists from Domains #1275

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

Merged
samuelpmishLLNL merged 13 commits into from
Nov 27, 2024
Merged

Fix acquisition of dof lists from Domains #1275

Show file tree
Hide file tree
Changes from 12 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
87fa6a8
Fix implementations of Enitre{Domain,Boundary} by removing repeated code
btalamini Nov 19, 2024
5a0f611
Put in test that catches bug in set operations, add tests for Entire*…
btalamini Nov 20, 2024
fff2744
Add a setter method to populate element lists to make object state st…
btalamini Nov 21, 2024
36c8669
Make setter method for adding boundary elements
btalamini Nov 21, 2024
8a9a11a
Simplify down to just one adder helper function
btalamini Nov 21, 2024
aee0ef3
Fix set operations
btalamini Nov 22, 2024
8506612
Refactor set operation application to avoid repetition
btalamini Nov 22, 2024
3052767
Make function to add lists of elements at once
btalamini Nov 22, 2024
d29960d
Refactor more repeated code
btalamini Nov 22, 2024
0e201b4
Add test of finding dofs on 3D elements
btalamini Nov 22, 2024
17f93be
Write docstrings
btalamini Nov 22, 2024
9cb9051
Make style
btalamini Nov 22, 2024
d4d79eb
fix bug when creating domains of boundary elements, delete vertex_ids…
samuelpmishLLNL Nov 26, 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
170 changes: 83 additions & 87 deletions src/serac/numerics/functional/domain.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -117,13 +117,11 @@ static Domain domain_of_edges(const mfem::Mesh& mesh, std::function<T> 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);
}
}
}
Expand Down Expand Up @@ -194,12 +192,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);
}
}

Expand Down Expand Up @@ -258,31 +254,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;
Expand All @@ -305,6 +297,39 @@ Domain Domain::ofElements(const mfem::Mesh& mesh, std::function<bool(std::vector
return domain_of_elems<3>(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<int>& geom_ids, const std::vector<int>& 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<int>::size_type i = 0; i < geom_ids.size(); ++i) {
addElement(geom_ids[i], elem_ids[i], element_geometry);
}
}

///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////

Expand Down Expand Up @@ -347,22 +372,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(edge_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(edge_id, f, geom);
}
quad_id++;
break;
Expand All @@ -385,6 +407,9 @@ Domain Domain::ofBoundaryElements(const mfem::Mesh& mesh, std::function<bool(std
return domain_of_boundary_elems<3>(mesh, func);
}

///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////

mfem::Array<int> Domain::dof_list(mfem::FiniteElementSpace* fes) const
{
std::set<int> dof_ids;
Expand Down Expand Up @@ -460,71 +485,32 @@ mfem::Array<int> 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<vec2>, int) { return true; });
break;
case 3:
return Domain::ofElements(mesh, [](std::vector<vec3>, 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<vec2>, int) { return true; });
break;
case 3:
return Domain::ofBoundaryElements(mesh, [](std::vector<vec3>, 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
Expand Down Expand Up @@ -553,22 +539,32 @@ Domain set_operation(set_op op, const Domain& a, const Domain& b)

Domain output{a.mesh_, a.dim_};

using Ids = std::vector<int>;
auto apply_set_op = [&op](const Ids& x, const Ids& y) { return set_operation(op, x, y); };

if (output.dim_ == 0) {
output.vertex_ids_ = set_operation(op, a.vertex_ids_, b.vertex_ids_);
output.vertex_ids_ = apply_set_op(a.vertex_ids_, b.vertex_ids_);
}

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;
Expand Down
52 changes: 45 additions & 7 deletions src/serac/numerics/functional/domain.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,14 +37,37 @@ 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
std::vector<int> vertex_ids_;

///@{
/// @name ElementIds
/// Indices of elements contained in the domain.
/// The first set, (vertex_ids_, edge_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.
///
/// @cond
std::vector<int> vertex_ids_;
/// 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.
std::vector<int> edge_ids_;
std::vector<int> tri_ids_;
std::vector<int> quad_ids_;
Expand All @@ -56,7 +79,7 @@ struct Domain {
std::vector<int> mfem_quad_ids_;
std::vector<int> mfem_tet_ids_;
std::vector<int> mfem_hex_ids_;
/// @endcond
///@}

Domain(const mfem::Mesh& m, int d, Type type = Domain::Type::Elements) : mesh_(m), dim_(d), type_(type) {}

Expand Down Expand Up @@ -132,6 +155,21 @@ 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.
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<int>& geom_id, const std::vector<int>& elem_id,
mfem::Geometry::Type element_geometry);
};

/// @brief constructs a domain from all the elements in a mesh
Expand Down
Loading