Skip to content

Commit

Permalink
Merge pull request #1275 from LLNL/talamini/bugfix/domain
Browse files Browse the repository at this point in the history 
Fix acquisition of dof lists from Domains
  • Loading branch information
@samuelpmishLLNL
samuelpmishLLNL authored Nov 27, 2024
2 parents 2932e9e + d4d79eb commit ad3c70a
Show file tree
Hide file tree
Showing 3 changed files with 254 additions and 193 deletions.
208 changes: 81 additions & 127 deletions src/serac/numerics/functional/domain.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,44 +48,6 @@ std::vector<tensor<double, d>> gather(const mfem::Vector& coordinates, mfem::Arr
return x;
}

template <int d>
static Domain domain_of_vertices(const mfem::Mesh& mesh, std::function<bool(tensor<double, d>)> 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<double, d> 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<bool(vec2)> func)
{
return domain_of_vertices(mesh, func);
}

Domain Domain::ofVertices(const mfem::Mesh& mesh, std::function<bool(vec3)> func)
{
return domain_of_vertices(mesh, func);
}

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

Expand Down Expand Up @@ -117,13 +79,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 +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);
}
}

Expand Down Expand Up @@ -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;
Expand All @@ -305,6 +259,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 +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;
Expand All @@ -385,6 +369,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 +447,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 +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<int>;
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;
Expand Down
51 changes: 45 additions & 6 deletions src/serac/numerics/functional/domain.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<int> vertex_ids_;
std::vector<int> edge_ids_;
std::vector<int> tri_ids_;
std::vector<int> quad_ids_;
Expand All @@ -58,6 +82,7 @@ struct Domain {
std::vector<int> 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) {}

/**
Expand Down Expand Up @@ -120,7 +145,6 @@ 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_;
Expand All @@ -132,6 +156,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

0 comments on commit ad3c70a

Please sign in to comment.