Add mesh generator that checks if a set of blocks is enclosed by a se…

…t of sidesets (#20621)

framework/doc/content/source/meshgenerators/SidesetsEncloseBlocks.md

@@ -0,0 +1,17 @@
+# SidesetsEncloseBlocks
+
+!syntax description /Mesh/SidesetsEncloseBlocks
+
+This mesh generator ensures that the blocks provided by the
+[!param](/meshgenerators/SidesetsEncloseBlocks/block)
+parameter are  enclosed by the sidesets provided by [!param](/meshgenerators/SidesetsEncloseBlocks/boundary).
+If that is not the case, the behavior of the mesh generator depends on whether
+[!param](/meshgenerators/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

framework/include/meshgenerators/SidesetsEncloseBlocks.h

@@ -0,0 +1,30 @@
+//* 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:
+  std::unique_ptr<MeshBase> & _input;
+};

framework/src/meshgenerators/SidesetsEncloseBlocks.C

@@ -0,0 +1,150 @@
+//* 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);
+
+defineLegacyParams(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);
+
+  // Get a reference to our BoundaryInfo object for later use
+  BoundaryInfo & boundary_info = mesh->get_boundary_info();
+
+  // get a list of all sides; vector of tuples (elem, loc_side, side_set)
+  auto side_list = boundary_info.build_active_side_list();
+
+  // 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);
+
+    // make sure that _sideset exists
+    if (new_sideset_id == BoundaryInfo::invalid_id)
+      paramError("new_boundary", "Not a valid boundary");
+  }
+
+  // get blocks
+  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 _sideset_name can be converted to integer it's interpreted
+  // as 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);
+
+    // make sure that _sideset exists
+    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)))
+    {
+      // loop through sides
+      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!!
+        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()));
+          // is intersection emtpy?
+          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);
+}

test/tests/meshgenerators/sideset_enclose_block/sideset_enclose_blocks.i

@@ -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
+  []
+[]

test/tests/meshgenerators/sideset_enclose_block/tests

@@ -0,0 +1,34 @@
+[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'
+    recover = false
+
+    [check_no_error]
+      type = 'Exodiff'
+      input = 'sideset_enclose_blocks.i'
+      exodiff = 'check_no_error.e'
+      detail = '.'
+      cli_args = '--mesh-only check_no_error.e'
+    []
+
+    [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.'
+    []
+
+    [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.'
+    []
+  []
+[]