Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
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
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
# SidesetsEncloseBlocks

!syntax description /Mesh/SidesetsEncloseBlocks

This mesh generator ensures that the blocks provided by the
[!param](/Mesh/SidesetsEncloseBlocks/block)
parameter are enclosed by the sidesets provided by [!param](/Mesh/SidesetsEncloseBlocks/boundary).
If that is not the case, the behavior of the mesh generator depends on whether
[!param](/Mesh/SidesetsEncloseBlocks/new_boundary) is provided.
If it is not provided, then an error is thrown. If it is provided, then the side
that is not covered is added to the new boundary.

!syntax parameters /Mesh/SidesetsEncloseBlocks

!syntax inputs /Mesh/SidesetsEncloseBlocks

!syntax children /Mesh/SidesetsEncloseBlocks
31 changes: 31 additions & 0 deletions framework/include/meshgenerators/SidesetsEncloseBlocks.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "MeshGenerator.h"

/**
* MeshGenerator that checks if a set of blocks is enclosed by a set of sidesets.
* It can either add sides that are not covered by a sideset by a new sidesets or
* error out.
*/
class SidesetsEncloseBlocks : public MeshGenerator
{
public:
static InputParameters validParams();

SidesetsEncloseBlocks(const InputParameters & parameters);

std::unique_ptr<MeshBase> generate() override;

protected:
/// the mesh that is passed from the meshgen executed before this meshgen
std::unique_ptr<MeshBase> & _input;
};
140 changes: 140 additions & 0 deletions framework/src/meshgenerators/SidesetsEncloseBlocks.C
Original file line number Diff line number Diff line change
@@ -0,0 +1,140 @@
//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "SidesetsEncloseBlocks.h"
#include "InputParameters.h"
#include "MooseTypes.h"
#include "MooseMeshUtils.h"
#include "CastUniquePointer.h"

#include "libmesh/remote_elem.h"

registerMooseObject("MooseApp", SidesetsEncloseBlocks);

InputParameters
SidesetsEncloseBlocks::validParams()
{
InputParameters params = MeshGenerator::validParams();

params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
params.addParam<std::vector<SubdomainName>>(
"block", "The set of blocks that are checked if they are enclosed by boundary");
params.addParam<std::vector<BoundaryName>>("boundary",
"The name of the boundaries enclosing the ");
params.addParam<BoundaryName>("new_boundary", "The name of the boundary to create");
params.addClassDescription(
"MeshGenerator that checks if a set of blocks is enclosed by a set of sidesets."
"It can either add sides that are not covered by a sideset by a new sidesets or"
"error out.");

return params;
}

SidesetsEncloseBlocks::SidesetsEncloseBlocks(const InputParameters & parameters)
: MeshGenerator(parameters), _input(getMesh("input"))
{
}

std::unique_ptr<MeshBase>
SidesetsEncloseBlocks::generate()
{
std::unique_ptr<MeshBase> mesh = std::move(_input);

BoundaryInfo & boundary_info = mesh->get_boundary_info();

// error on finding a side that is not covered
bool error_out = !isParamValid("new_boundary");
boundary_id_type new_sideset_id;
if (!error_out)
{
BoundaryName new_boundary = getParam<BoundaryName>("new_boundary");
std::stringstream ss;
ss << new_boundary;
ss >> new_sideset_id;
if (ss.fail())
new_sideset_id = boundary_info.get_id_by_name(new_boundary);

if (new_sideset_id == BoundaryInfo::invalid_id)
paramError("new_boundary", "Not a valid boundary");
}

std::vector<subdomain_id_type> vec_block_ids =
MooseMeshUtils::getSubdomainIDs(*mesh, getParam<std::vector<SubdomainName>>("block"));
std::set<subdomain_id_type> blk_ids(vec_block_ids.begin(), vec_block_ids.end());

// get boundaries
// check if the provided sideset name is actually a sideset id
// if it can be converted to integer it's interpreted
// as a sideset id
auto boundary_name_vec = getParam<std::vector<BoundaryName>>("boundary");
std::vector<boundary_id_type> bnd_ids_vec;
bnd_ids_vec.reserve(boundary_name_vec.size());
std::set<boundary_id_type> bnd_ids;
for (auto & bnd_name : boundary_name_vec)
{
std::stringstream ss;
boundary_id_type sideset_id;
ss << bnd_name;
ss >> sideset_id;
if (ss.fail())
sideset_id = boundary_info.get_id_by_name(bnd_name);

if (sideset_id == BoundaryInfo::invalid_id)
paramError("boundary", "Not a valid boundary");
bnd_ids_vec.push_back(sideset_id);
bnd_ids.insert(sideset_id);
}

// loop over all elements in blocks and for each elem over each side
// and check the neighbors
for (auto & block_id : blk_ids)
{
for (const Elem * elem : as_range(mesh->active_local_subdomain_elements_begin(block_id),
mesh->active_local_subdomain_elements_end(block_id)))
{
Comment on lines +96 to +100
Copy link
Member

Choose a reason for hiding this comment

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

I think I would do something like the following, but I would have @roystgnr confirm

Suggested change
for (auto & block_id : blk_ids)
{
for (const Elem * elem : as_range(mesh->active_local_subdomain_elements_begin(block_id),
mesh->active_local_subdomain_elements_end(block_id)))
{
for (const Elem * elem : as_range(mesh->local_level_elements_begin(0),
mesh->local_level_elements_end(0)))
{
if (!blk_ids.count(elem->subdomain_id()))
continue;

I would only bother looping over level 0 because in fact we error if you attempt to add sideset info for a non-level 0 element

void BoundaryInfo::add_side(const Elem * elem,
                            const unsigned short int side,
                            const boundary_id_type id)
{
  libmesh_assert(elem);

  // Only add BCs for level-0 elements.
  libmesh_assert_equal_to (elem->level(), 0);

Copy link
Contributor

Choose a reason for hiding this comment

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

"You can't add sideset info for a non-level-0 element" is more of a missing feature than an API design. We're hoping to fix that in the near future, once @fdkong or I have time to look at libMesh/libmesh#3094 again.

"If you add boundary info to a parent element, that automatically associates all its descendant elements sharing that side with that boundary id", on the other hand, is an API design we'll be leaving unchanged.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

So say if I hit an element that was refined by AMR and I try to add one of its sides to a sideset, and I run in devel or dbg, then I would hit is assert ( libmesh_assert_equal_to (elem->level(), 0);)?

In that case, the code cannot safely work for AMR and potentially not for distributed mesh either. I know how to check for distributed mesh and error out. For AMR, can I check if the active element is the ultimate parent or something like that?

Copy link
Contributor

Choose a reason for hiding this comment

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

Hmm... checking if an active element is level-0 (the ultimate parent) isn't a good AMR check (a uniformly refined mesh will have active elements at higher levels), but it might be a good check for your situation here, for now: when you're looping over active elements just error out if you're not at level 0. Even on that hypothetical uniformly refined mesh, you won't currently be allowed to add boundary ids to child sides.

If we didn't have libMesh/libmesh#3094 waiting in the wings I'd say it would be worth a little effort to try to handle refined meshes where child elements do all belong to their parents' subdomains, but instead the best thing for you to do might be to just handle the unrefined mesh case for now but tag yourself on that PR so we don't forget to notify you when handling the refined case too becomes an option.

for (unsigned int j = 0; j < elem->n_sides(); ++j)
{
const Elem * neigh = elem->neighbor_ptr(j);

// is this an outside boundary to blocks?
// NOTE: the next line is NOT sufficient for AMR!!
Copy link
Contributor Author

Choose a reason for hiding this comment

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

can I have some help from maybe @lindsayad. I don't think the check I am doing will work if we have hanging nodes. How do I generalize this correctly?

Copy link
Contributor

Choose a reason for hiding this comment

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

Are you trying to support the case where you have subdomains that differ on different children of the same parent element? libMesh supports that case, but unfortunately since we don't yet support sidesets that aren't on sides of parent elements, you won't yet be able to support that case here. For now it probably makes sense to just worry about level-0.

I'd worry more about distributed mesh support, personally. For the sides in between ghost elements and remote elements on any particular processor's part of the mesh, you don't have enough information to tell whether you need a new boundary id there or not without communication. The solutions to this are "exit with an error if the mesh is distributed", "throw in a MeshSerializer to temporarily un-distribute it" or "write that communication routine", ranked from least to most by both difficulty and usefulness.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

My concern was this case

Screen Shot 2022-03-30 at 12 29 43 PM

Both elements are in the same block. But between elem 2 is refined. The side is checked
with this line of code:

bool is_outer_bnd = !neigh || blk_ids.find(neigh->subdomain_id()) == blk_ids.end();

My fear was that somehow neigh == nullptr in that case. I now think that neigh is not
nullptr but that it is simply not an active element. Is that correct @roystgnr?

I think that is different than what you thought I
wanted. Is that the case you described:
Screen Shot 2022-03-30 at 12 35 23 PM
i.e., a single parent element has children with different blocks. That is too exotic of a case for me to care about.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Actually if I think about it. What subdomain is returned for this case:
Screen Shot 2022-03-30 at 12 38 08 PM
I am in elem 1 and for the indicated side I get the neighbor pointer. Then I get the subdomain id. What is that subdomain id? (B1, B2, B3)?

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'd worry more about distributed mesh support, personally. For the sides in between ghost elements and remote elements on any particular processor's part of the mesh, you don't have enough information to tell whether you need a new boundary id there or not without communication. The solutions to this are "exit with an error if the mesh is distributed", "throw in a MeshSerializer to temporarily un-distribute it" or "write that communication routine", ranked from least to most by both difficulty and usefulness.

I am not sure I understand why that is. I don't think I ever access sides of ghosted elements because my loops are only over active local elems.

Copy link
Contributor

Choose a reason for hiding this comment

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

In the first case, where I'll assume you got four children by refining Elem 2 and you haven't edited any of their subdomain ids since, I think you're safe: Elem 2 still has the same subdomain as all it's children, and Elem 2 is still the neighbor of Elem 1, it's just that elem2->active() would now be false for it.

In the second case? You're not at all safe: for all you know, you might get B4 for the subdomain id! Elem::subdomain_id() just returns the subdomain this element, regardless of whether it's active or not, and regardless of what ids its descendants might have. If, for instance, elem2 is on block 4, and you refine it and set its children to live on blocks 1,2, and 3, then elem2 is still on block 4! You can get its set of active descendants (or its set of active descendants on the side neighboring elem1; that's a natural enough task that we have an API for it) and check each of their subdomain_id() values individually, but we won't do anything to try to keep the parent in sync here; in cases like this where you give different subdomains to different children it's impossible to keep the parent in sync with all of them.

bool is_outer_bnd = !neigh || blk_ids.find(neigh->subdomain_id()) == blk_ids.end();
if (is_outer_bnd)
{
// get all boundary ids of this side, then compare the set of these boundary_ids
// to the set of boundary_ids provided to the mesh generator; if the intersection
// is empty then this side is NOT convered
std::vector<boundary_id_type> side_boundary_ids_vec;
boundary_info.raw_boundary_ids(elem, j, side_boundary_ids_vec);

std::set<boundary_id_type> intersection;
std::set_intersection(side_boundary_ids_vec.begin(),
side_boundary_ids_vec.end(),
bnd_ids_vec.begin(),
bnd_ids_vec.end(),
std::inserter(intersection, intersection.end()));

if (intersection.size() == 0)
{
if (error_out)
mooseError("Element id ",
elem->id(),
" side ",
j,
" is external and not covered by specified boundaries.");
else
boundary_info.add_side(elem, j, new_sideset_id);
}
}
}
}
}

return dynamic_pointer_cast<MeshBase>(mesh);
}
Binary file not shown.
Binary file not shown.
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
[Mesh]
[cmg]
type = CartesianMeshGenerator
dim = 2
dx = '1 1 1 1'
dy = '1 1 1 1'
subdomain_id = '1 2 1 2
1 3 4 4
2 5 10 1
1 2 2 2'
[]

[add_sides]
type = SideSetsBetweenSubdomainsGenerator
primary_block = '3 4 5 10'
paired_block = '1 2'
new_boundary = 'interior'
input = cmg
[]

[enclosed]
type = SidesetsEncloseBlocks
block = '3 4 5 10'
boundary = '1 4'
input = add_sides
[]
[]
35 changes: 35 additions & 0 deletions test/tests/meshgenerators/sideset_enclose_block/tests
Original file line number Diff line number Diff line change
@@ -0,0 +1,35 @@
[Tests]
design = 'meshgenerators/SidesetsEncloseBlocks.md'
issues = '#20621'

[sidesets_enclosed_blocks]
requirement = 'The system shall be able to check of a set of blocks is enclosed by a set of sidesets'
[check_no_error]
type = 'Exodiff'
input = 'sideset_enclose_blocks.i'
exodiff = 'check_no_error.e'
detail = '.'
cli_args = '--mesh-only check_no_error.e'
recover = false
[]

[check_error]
type = RunException
input = 'sideset_enclose_blocks.i'
expect_err = "Element id 5 side 0 is external and not covered by specified boundaries."
cli_args = 'Mesh/add_sides/paired_block="1" --mesh-only'
detail = ' and error out if this condition is not met.'
recover = false
[]

[add_new_boundary]
type = 'Exodiff'
input = 'sideset_enclose_blocks.i'
exodiff = 'add_new_boundary.e'
max_parallel = 1
cli_args = 'Mesh/add_sides/paired_block="1" Mesh/enclosed/new_boundary=12 --mesh-only add_new_boundary.e'
detail = ' and upon user request add missing sides into a specified input.'
recover = false
[]
[]
[]