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

Allow boundary be associated with child elements #3094

Open
wants to merge 3 commits into
base: devel
Choose a base branch
from
Open

Allow boundary be associated with child elements #3094

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
82491b1
Allow boundary be associated with child elements
fdkong Nov 29, 2021
471836f
Added APIs to allow users set the flag: _children_on_boundary
fdkong Feb 23, 2022
4be946d
Made BoundaryInfo::boundary_ids for _children_on_boundary
fdkong Feb 24, 2022
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
19 changes: 19 additions & 0 deletions include/mesh/boundary_info.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -874,6 +874,18 @@ class BoundaryInfo : public ParallelObject
const std::multimap<const Elem *, std::pair<unsigned short int, boundary_id_type>> & get_sideset_map() const
{ return _boundary_side_id; }

/**
* \returns Whether or not there are some children on boundary sides
*/
bool is_children_on_boundary_side() const
{ return _children_on_boundary; }

/**
* Whether or not to allow set boundary sides on children elements
*/
void allow_children_on_boundary_side(const bool children_on_boundary)
{ _children_on_boundary = children_on_boundary; }

private:

/**
Expand Down Expand Up @@ -924,6 +936,13 @@ class BoundaryInfo : public ParallelObject
std::pair<unsigned short int, boundary_id_type>>
_boundary_side_id;

/*
* Whether or not children elements are associated to any boundary
* It is false by default. The flag will be turnned on if add_side
* function is called with a child element
*/
bool _children_on_boundary;

/**
* A collection of user-specified boundary ids for sides, edges, nodes,
* and shell faces.
Expand Down
85 changes: 63 additions & 22 deletions src/mesh/boundary_info.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,8 @@ const boundary_id_type BoundaryInfo::invalid_id = -123;
// BoundaryInfo functions
BoundaryInfo::BoundaryInfo(MeshBase & m) :
ParallelObject(m.comm()),
_mesh (&m)
_mesh (&m),
_children_on_boundary(false)
{
}

Expand Down Expand Up @@ -950,8 +951,11 @@ void BoundaryInfo::add_side(const Elem * elem,
{
libmesh_assert(elem);

// Only add BCs for level-0 elements.
libmesh_assert_equal_to (elem->level(), 0);
// Users try to mark boundary on child elements
// If this happens, we will allow users to remove
// side from child elements as well
if (elem->level())
_children_on_boundary = true;

libmesh_error_msg_if(id == invalid_id, "ERROR: You may not set a boundary ID of "
<< invalid_id
Expand Down Expand Up @@ -979,8 +983,11 @@ void BoundaryInfo::add_side(const Elem * elem,

libmesh_assert(elem);

// Only add BCs for level-0 elements.
libmesh_assert_equal_to (elem->level(), 0);
// Users try to mark boundary on child elements
// If this happens, we will allow users to remove
// side from child elements as well
if (elem->level())
_children_on_boundary = true;

// Don't add the same ID twice
auto bounds = _boundary_side_id.equal_range(elem);
Expand Down Expand Up @@ -1228,8 +1235,15 @@ void BoundaryInfo::boundary_ids (const Elem * const elem,
// Clear out any previous contents
vec_to_fill.clear();

// Only level-0 elements store BCs. If this is not a level-0
// element get its level-0 parent and infer the BCs.
// Search BC on the current element. If we find anything, we should return
for (const auto & pr : as_range(_boundary_side_id.equal_range(elem)))
if (pr.second.first == side)
vec_to_fill.push_back(pr.second.second);

if (vec_to_fill.size())
return;

// We check the top parent now
const Elem * searched_elem = elem;
if (elem->level() != 0)
{
Expand Down Expand Up @@ -1276,7 +1290,7 @@ void BoundaryInfo::raw_boundary_ids (const Elem * const elem,
vec_to_fill.clear();

// Only level-0 elements store BCs.
if (elem->parent())
if (elem->parent() && !_children_on_boundary)
return;

// Check each element in the range to see if its side matches the requested side.
Expand Down Expand Up @@ -1426,8 +1440,8 @@ void BoundaryInfo::remove_side (const Elem * elem,
{
libmesh_assert(elem);

// Only level 0 elements are stored in BoundaryInfo.
libmesh_assert_equal_to (elem->level(), 0);
// Only level 0 elements unless the flag "_children_on_boundary" is on.
libmesh_assert(elem->level() == 0 || _children_on_boundary);

// Erase (elem, side, *) entries from map.
erase_if(_boundary_side_id, elem,
Expand All @@ -1443,6 +1457,9 @@ void BoundaryInfo::remove_side (const Elem * elem,
{
libmesh_assert(elem);

// Only level 0 elements unless the flag "_children_on_boundary" is on.
libmesh_assert(elem->level() == 0 || _children_on_boundary);

// Erase (elem, side, id) entries from map.
erase_if(_boundary_side_id, elem,
[side, id](decltype(_boundary_side_id)::mapped_type & pr)
Expand Down Expand Up @@ -1489,12 +1506,23 @@ void BoundaryInfo::remove_id (boundary_id_type id)
unsigned int BoundaryInfo::side_with_boundary_id(const Elem * const elem,
const boundary_id_type boundary_id_in) const
{
const Elem * searched_elem = elem;
std::vector<const Elem *> searched_elem_vec;
// If elem has boundary information, we return that when
// the flag "_children_on_boundary" is on
if (_children_on_boundary)
searched_elem_vec.push_back(elem);
// Otherwise, we return boundary information of its
// parent if any
if (elem->level() != 0)
searched_elem = elem->top_parent();
searched_elem_vec.push_back(elem->top_parent());
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's not enough to just test elem and elem->top_parent(), is it? Suppose we add a boundary side to a child element, then refine that element. The grandchildren elements on that side should still have that boundary id, shouldn't they? We'll want to loop upward parent by parent in the _chidren_on_boundary case.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For my use case, it is enough. Top parent is assigned a boundary when reading in a mesh. Once the boundary is moving, boundaries will be on active elements. We literally delete everything else.

Should we allow users to leave boundary everywhere? I can make that change

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I feel like if the current element has a boundary, and then parents should not have the same boundary ID. Otherwise which one we should choose for users?

So do we need to do some checks?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One person's "Works for me" is another person's "Booby trap". Refining an element gives you children who have the same boundary ids on the sides they share with their parent; for that to suddenly fail is not a bug libMesh currently has and we're not about to introduce it now.

I feel like if the current element has a boundary, and then parents should not have the same boundary ID. Otherwise which one we should choose for users?

If they're all the same then there's no choice; we return that id. If they differ, then we return true if we match any of them (in the has_boundary_id test case) or we add them all to the vector we're filling (in the boundary_ids test case). Say element p has boundary id b1 on a side. We refine it, and on the shared side child c now also correctly has boundary id b1. What should happen when we add b2 to that side of c? That shouldn't imply removing b1. When querying for a complete set of boundary ids we'll need to remember to query p (and on up the tree) even if we find one or more ids assigned to lower levels.

So do we need to do some checks?

I think the one place we're missing here where it gets a little interesting is coarsening.

In general, when we coarsen an element, the boundary data attached to the now-vanishing children needs to disappear, not to become dangling pointers.

When we try to coarsen the mesh down to a parent element which has had the same boundary id added to that side of all its children. If the id hasn't been explicitly added to the parent, it won't be on the parent. But if that side was entirely on the boundary before coarsening it should be on entirely on the boundary after! In that case we should be adding the id to the parent as we remove it from the children.

And the only place that gets really interesting: what happens when some of the children have an id and others do not? I could see an argument for "the parent then gets the id from any of its children" or for "the parent doesn't get the id" or for "the parent gets the id if a majority (rounding up? down?) of its children have it", and that's all fine because when we're coarsening we expect our data to be getting coarser, but whatever we do we ought to document the hell out of it and add a few unit tests.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, wait, I could see one more interesting question: if we do want to enforce "a parent side gets ids that were given to one/most of the children sides it contains" (we may) or "gets ids that were given to all of the children sides" (we do?), do we want to add that id to the parent preemptively, when it's added to (enough/all of?) the children, before coarsening? That might be useful for users to be able to query directly without ever having to manually set ids on non-active elements. But then on the other hand if we do that then we need to do the same for removal of ids, for consistency, and that very consistency would a source of inflexibility for users who might prefer a different behavior than whatever we choose.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @roystgnr , thanks for your review! I am trying to take over the PR from @fdkong so that this effort does not go down in the drain. However I am super confused... What is the agreement that you guys arrive at? - Revert the changes here and add a documentation about not returning descendants' boundary ids from this function? Thanks :)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

IIRC the code functionalities we absolutely need to add are "search all ancestors when filling out children's boundary ids" and "transfer an id up to an ancestor when coarsening away children that all have that id on that ancestor's side". The bit we need to decide upon is what to do when coarsening away children only some of whom have an id on a particular side, and for that I had three suggestions, but I'd even be content with the laziest of the three, "the id gets lost there", so long as it's documented well.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @roystgnr ! Please correct me if I am wrong:

search all ancestors when filling out children's boundary ids

so this is to avoid the case of not returning the element side if any of its "intermediate" ancestors have the boundary ID while its top parent does not?

transfer an id up to an ancestor when coarsening away children that all have that id on that ancestor's side

guess this should be done in the function BoundaryInfo::add_sides, not in this function?

what to do when coarsening away children only some of whom have an id on a particular side, and for that I had three suggestions, but I'd even be content with the laziest of the three, "the id gets lost there", so long as it's documented well.

not very sure where and how to keep track of (or throw an exception/warning..?) about children elements of different IDs being coarsened away... would you mind pointing me to the function where this should be added?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

to avoid the case of not returning the element side if any of its "intermediate" ancestors have the boundary ID while its top parent does not?

Exactly.

guess this should be done in the function BoundaryInfo::add_sides, not in this function?

Hmm... actually, I'd say it's not safe to do it until the actual coarsening happens. If we try to do it any earlier, then something as simple as adding a boundary id and removing it again could see the removal fail, because that addition was the last child sharing a parent side and we "helpfully" added the id to the parent. Or ... we could also propagate removals up to ancestors? That would make sense too: when the last child sharing a side has an id added there, we add it to the parent side; when any child sharing a side has an id removed we remove it from the parent (and so recursively from any other ancestors that have it)? That gets really tricky, though. The desired behavior is easy enough to describe, but the implementation details...

not very sure where and how to keep track of (or throw an exception/warning..?) about children elements of different IDs being coarsened away

We don't throw warnings or exceptions about losing information when coarsening; that's kind of implied by the concept of coarsening. It just needs to be in the BoundaryInfo documentation somewhere.

would you mind pointing me to the function where this should be added?

MeshRefinement::coarsen_elements() could handle transferring boundary ids to newly-reactivated coarse parents, maybe?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Cool! Really appreciate your suggestions here! Let me look into this a bit more and I may (will) bother you again with some questions about the implementation details.

else if (!_children_on_boundary)
searched_elem_vec.push_back(elem);

// elem may have zero or multiple occurrences
for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
for (auto it = searched_elem_vec.begin(); it != searched_elem_vec.end(); ++it)
{
const Elem * searched_elem = *it;
// elem may have zero or multiple occurrences
for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
{
// if this is true we found the requested boundary_id
// of the element and want to return the side
Expand Down Expand Up @@ -1525,7 +1553,8 @@ unsigned int BoundaryInfo::side_with_boundary_id(const Elem * const elem,
if (!p)
return side;
}
}
}
}

// if we get here, we found elem in the data structure but not
// the requested boundary id, so return the default value
Expand All @@ -1539,12 +1568,22 @@ BoundaryInfo::sides_with_boundary_id(const Elem * const elem,
{
std::vector<unsigned int> returnval;

const Elem * searched_elem = elem;
std::vector<const Elem *> searched_elem_vec;
// If elem has boundary information, that is part of return when
// the flag "_children_on_boundary" is on
if (_children_on_boundary)
searched_elem_vec.push_back(elem);
// Return boundary information of its parent as well
if (elem->level() != 0)
searched_elem = elem->top_parent();
searched_elem_vec.push_back(elem->top_parent());
else if (!_children_on_boundary)
searched_elem_vec.push_back(elem);

// elem may have zero or multiple occurrences
for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
for (auto it = searched_elem_vec.begin(); it != searched_elem_vec.end(); ++it)
{
const Elem * searched_elem = *it;
// elem may have zero or multiple occurrences
for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
{
// if this is true we found the requested boundary_id
// of the element and want to return the side
Expand Down Expand Up @@ -1579,7 +1618,7 @@ BoundaryInfo::sides_with_boundary_id(const Elem * const elem,
returnval.push_back(side);
}
}

}
return returnval;
}

Expand Down Expand Up @@ -1794,7 +1833,8 @@ BoundaryInfo::build_node_list_from_side_list()
// Need to loop over the sides of any possible children
std::vector<const Elem *> family;
#ifdef LIBMESH_ENABLE_AMR
elem->active_family_tree_by_side (family, id_pair.first);
if (!elem->subactive())
elem->active_family_tree_by_side (family, id_pair.first);
#else
family.push_back(elem);
#endif
Expand Down Expand Up @@ -2156,7 +2196,8 @@ BoundaryInfo::build_active_side_list () const
// Loop over the sides of possible children
std::vector<const Elem *> family;
#ifdef LIBMESH_ENABLE_AMR
elem->active_family_tree_by_side(family, id_pair.first);
if (!elem->subactive())
elem->active_family_tree_by_side(family, id_pair.first);
#else
family.push_back(elem);
#endif
Expand Down