Add a mesh generator for creating a polyline that follows a nodeset

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

@@ -0,0 +1,28 @@
+# PolyLineMeshFollowingNodeSetGenerator
+
+!syntax description /Mesh/PolyLineMeshFollowingNodeSetGenerator
+
+Using the `PolyLineMeshFollowingNodeSetGenerator` object from within the
+[Mesh](/Mesh/index.md) block of the input file will construct an open
+or closed (looped) one-dimensional manifold of Edge elements.
+
+The polyline is generated iteratively, starting from the [!param](/Mesh/PolyLineMeshFollowingNodeSetGenerator/starting_point),
+taking the first step following the [!param](/Mesh/PolyLineMeshFollowingNodeSetGenerator/starting_direction):
+The following heuristics are implemented for following the nodeset:
+
+- Sphere-based nodeset centroid search
+    - the current point is moved by [!param](/Mesh/PolyLineMeshFollowingNodeSetGenerator/dx) towards the current direction
+    - the local centroid of all nodes in the target nodeset and closer than [!param](/Mesh/PolyLineMeshFollowingNodeSetGenerator/search_radius) to the current point is found
+    - the current direction is updated to be from the current point to the centroid
+    - the displacement of the current point is annulled, and replaced by a displacement of magnitude [!param](/Mesh/PolyLineMeshFollowingNodeSetGenerator/dx) towards the centroid
+- Ignoring nodes behind, activated with the [!param](/Mesh/PolyLineMeshFollowingNodeSetGenerator/ignore_nodes_behind) parameter
+    - only nodes from the nodeset that are also ahead of the current point with the current direction are considered
+
+!alert note
+Because of the use of heuristics rather than a global solve, the 1D polyline solution is highly dependent on the parameters chosen.
+The user should tune these parameters to fit their needs. Notably, the number of edges can be reduced to prevent the polyline from
+backtracking after having reached the end of the nodeset.
+
+!syntax parameters /Mesh/PolyLineMeshFollowingNodeSetGenerator
+
+!syntax inputs /Mesh/PolyLineMeshFollowingNodeSetGenerator

framework/include/meshgenerators/PolyLineMeshFollowingNodeSetGenerator.h

@@ -0,0 +1,52 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* 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"
+
+/**
+ * Generates a polyline (open ended or looped) of Edge elements
+ * by marching along a nodeeset and trying to be as close as possible to
+ * the nodes of the nodeset
+ */
+class PolyLineMeshFollowingNodeSetGenerator : public MeshGenerator
+{
+public:
+  static InputParameters validParams();
+
+  PolyLineMeshFollowingNodeSetGenerator(const InputParameters & parameters);
+
+  std::unique_ptr<MeshBase> generate() override;
+
+protected:
+  // Input mesh with the nodeset
+  std::unique_ptr<MeshBase> & _input;
+
+  /// Starting point of the polyline
+  const Point _starting_point;
+  /// Starting direction for the polyline
+  const Point _starting_direction;
+  /// Whether to only look in front of the polyline for nodes in the nodeset
+  const bool _ignore_nodes_behind;
+  /// Whether edges should form a closed loop. Will error if the nodeset does not loop back on itself
+  const bool _loop;
+
+  /// Subdomain name to assign to the polyline edge elements
+  const SubdomainName _line_subdomain;
+  /// Boundary names to assign to (non-looped) polyline start and end
+  const BoundaryName _start_boundary, _end_boundary;
+
+  /// Approximate spacing between nodes
+  const Real _dx;
+  /// How many Edge elements to build between each point pair
+  const unsigned int _num_edges_between_points;
+  /// Whether to output to console the mesh generation process
+  const bool _verbose;
+};

framework/src/meshgenerators/PolyLineMeshFollowingNodeSetGenerator.C

@@ -0,0 +1,211 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* 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 "PolyLineMeshFollowingNodeSetGenerator.h"
+
+#include "CastUniquePointer.h"
+#include "MooseMeshUtils.h"
+#include "MooseUtils.h"
+
+#include "libmesh/elem.h"
+#include "libmesh/int_range.h"
+#include "libmesh/unstructured_mesh.h"
+
+registerMooseObject("MooseApp", PolyLineMeshFollowingNodeSetGenerator);
+
+InputParameters
+PolyLineMeshFollowingNodeSetGenerator::validParams()
+{
+  InputParameters params = MeshGenerator::validParams();
+
+  // Path parameters
+  params.addRequiredParam<Point>("starting_point", "Starting point for the polyline");
+  params.addRequiredParam<Point>("starting_direction",
+                                 "Starting value for the direction of the line");
+  params.addRequiredParam<MeshGeneratorName>("input", "The mesh we get the sideset from");
+  params.addRequiredParam<BoundaryName>("nodeset", "Nodeset to follow to form the polyline");
+  params.addRequiredParam<Real>("search_radius",
+                                "Radius of the sphere used to find points in the nodeset");
+  params.addParam<bool>(
+      "ignore_nodes_behind",
+      false,
+      "Ignore nodes in the nodeset that are behind the current point in the polyline");
+  params.addParam<bool>("loop", false, "Whether edges should form a closed loop");
+
+  // Discretization parameters
+  // NOTE: we could have another dx as a path search parameter, and decouple the two options
+  params.addRequiredRangeCheckedParam<Real>("dx", "dx>0", "Approximate size of the edge elements");
+  params.addParam<unsigned int>(
+      "max_edges", 1000, "Maximum number of edges. Serves as a stopping criterion");
+  params.addParam<unsigned int>(
+      "num_edges_between_points", 1, "How many Edge elements to build between each point pair");
+
+  // Naming of result parameters
+  params.addParam<SubdomainName>("line_subdomain", "line", "Subdomain name for the line");
+  params.addParam<BoundaryName>(
+      "start_boundary", "start", "Boundary to assign to (non-looped) polyline start");
+  params.addParam<BoundaryName>(
+      "end_boundary", "end", "Boundary to assign to (non-looped) polyline end");
+
+  params.addParam<bool>(
+      "verbose",
+      false,
+      "whether to output additional information to console during the line generation");
+
+  params.addClassDescription(
+      "Generates a polyline (open ended or looped) of Edge elements by marching along a nodeeset "
+      "and trying to be as close as possible to the nodes of the nodeset");
+
+  return params;
+}
+
+PolyLineMeshFollowingNodeSetGenerator::PolyLineMeshFollowingNodeSetGenerator(
+    const InputParameters & parameters)
+  : MeshGenerator(parameters),
+    _input(getMesh("input")),
+    _starting_point(getParam<Point>("starting_point")),
+    _starting_direction(getParam<Point>("starting_direction")),
+    _ignore_nodes_behind(getParam<bool>("ignore_nodes_behind")),
+    _loop(getParam<bool>("loop")),
+    _line_subdomain(getParam<SubdomainName>("line_subdomain")),
+    _start_boundary(getParam<BoundaryName>("start_boundary")),
+    _end_boundary(getParam<BoundaryName>("end_boundary")),
+    _dx(getParam<Real>("dx")),
+    _num_edges_between_points(getParam<unsigned int>("num_edges_between_points")),
+    _verbose(getParam<bool>("verbose"))
+{
+  if (_loop && (isParamSetByUser("start_boundary") || isParamSetByUser("end_boundary")))
+    paramError("loop",
+               "Loop does not have a start or end boundary. These parameters must not be passed.");
+}
+
+std::unique_ptr<MeshBase>
+PolyLineMeshFollowingNodeSetGenerator::generate()
+{
+  auto uptr_mesh = buildMeshBaseObject();
+  MeshBase & mesh = *uptr_mesh;
+  std::unique_ptr<MeshBase> base_mesh = std::move(_input);
+  if (!base_mesh->is_serial())
+    paramError("input", "Input mesh must not be distributed");
+
+  // Get nodeset ID in input mesh
+  const auto nodeset_id =
+      MooseMeshUtils::getBoundaryID(getParam<BoundaryName>("nodeset"), *base_mesh);
+
+  const auto search_radius_sq = std::pow(getParam<Real>("search_radius"), 2);
+  const auto n_points = getParam<unsigned int>("max_edges");
+  Point current_point = _starting_point;
+  Point previous_direction = _starting_direction;
+  mesh.add_point(_starting_point, 0);
+
+  // Pre-find all points in the nodeset to follow
+  // TODO: Build a KNN tree to speed up the search later on
+  const auto all_nodeset_tuples = base_mesh->get_boundary_info().build_node_list(
+      libMesh::BoundaryInfo::NodeBCTupleSortBy::BOUNDARY_ID);
+  std::vector<dof_id_type> nodeset_nodes;
+  nodeset_nodes.reserve(all_nodeset_tuples.size() /
+                        base_mesh->get_boundary_info().n_boundary_ids());
+  for (const auto & tup : all_nodeset_tuples)
+    if (BoundaryID(std::get<1>(tup)) == nodeset_id)
+      nodeset_nodes.push_back(std::get<0>(tup));
+  if (_verbose)
+    _console << "Total number of nodes in nodeset " << getParam<BoundaryName>("nodeset") << ": "
+             << nodeset_nodes.size() << std::endl;
+
+  unsigned int n_segments = 0;
+  for (auto i : make_range(n_points))
+  {
+    // Move the point forward in the search direction
+    const auto previous_point = current_point;
+    current_point += previous_direction * _dx;
+
+    // Draw a sphere to find the next point as the barycenter of the nodes from the nodeset inside
+    // the sphere.
+    // NOTE: there are many heuristics we could use here. We could try to fit a cylinder
+    // to the nodeset if we know the nodeset represents a cylinder for example.
+    Point barycenter(0);
+    unsigned int n_sum = 0;
+    for (const auto n_id : nodeset_nodes)
+      if ((current_point - base_mesh->node_ref(n_id)).norm_sq() < search_radius_sq)
+      {
+        if (!_ignore_nodes_behind ||
+            ((base_mesh->node_ref(n_id) - current_point) * previous_direction >= 0))
+        {
+          barycenter += base_mesh->node_ref(n_id);
+          n_sum++;
+        }
+      }
+    if (n_sum > 0)
+      barycenter /= n_sum;
+    else if (!_ignore_nodes_behind)
+      mooseError("Did not find any nodes in the nodeset near the current point at: ",
+                 current_point);
+    else
+      barycenter = previous_point;
+
+    if (MooseUtils::absoluteFuzzyEqual((barycenter - previous_point).norm_sq(), 0))
+    {
+      mooseInfo("Barycenter did not move from ", barycenter, ". Returning!");
+      goto done_drawing;
+    }
+
+    // Compute new direction
+    const auto new_direction = (barycenter - previous_point).unit();
+    previous_direction = new_direction;
+
+    // Set the new point towards the barycenter
+    n_segments++;
+    // Note: this dx could be different than the dx used to search the barycenter
+    current_point = previous_point + _dx * new_direction;
+    mesh.add_point(current_point, (i + 1) * _num_edges_between_points);
+    if (_verbose)
+      _console << i << ": new point: " << current_point << " new direction " << previous_direction
+               << std::endl;
+
+    // Add the additional edges in between if requested
+    if (_num_edges_between_points > 1)
+    {
+      auto p = previous_point;
+      const Point pvec = (current_point - previous_point) / _num_edges_between_points;
+      for (auto j : make_range(1u, _num_edges_between_points))
+      {
+        p += pvec;
+        mesh.add_point(p, i * _num_edges_between_points + j);
+      }
+    }
+  }
+
+done_drawing:
+
+  const auto n_elem = n_segments * _num_edges_between_points - 1;
+  const auto max_nodes = n_segments * _num_edges_between_points - 1;
+  for (auto i : make_range(n_elem + _loop))
+  {
+    const auto ip1 = _loop ? (i + 1) % max_nodes : (i + 1);
+    auto elem = Elem::build(EDGE2);
+    elem->set_node(0, mesh.node_ptr(i));
+    elem->set_node(1, mesh.node_ptr(ip1));
+    elem->set_id() = i;
+    mesh.add_elem(std::move(elem));
+  }
+
+  // Add the starting and end boundary
+  if (!_loop)
+  {
+    BoundaryInfo & bi = mesh.get_boundary_info();
+    std::vector<BoundaryName> bdy_names{_start_boundary, _end_boundary};
+    std::vector<boundary_id_type> ids = MooseMeshUtils::getBoundaryIDs(mesh, bdy_names, true);
+    bi.add_side(mesh.elem_ptr(0), 0, ids[0]);
+    bi.add_side(mesh.elem_ptr(n_elem - 1), 1, ids[1]);
+  }
+
+  mesh.prepare_for_use();
+
+  return uptr_mesh;
+}

test/tests/meshgenerators/polyline_following_nodeset_mesh_generator/2d_channel.i

@@ -0,0 +1,50 @@
+[Mesh]
+  [gmg]
+    type = CartesianMeshGenerator
+    dim = 2
+    ix = '1 2 2 1'
+    dx = '1 1 1 1'
+    # channel has details in x direction, y-discretization helps guide the tube. This was
+    # not necessary for the cylindrical channel in 3D where the bends did not influence as much
+    # the direction of the polyline
+    iy = '5 4 4 2'
+    dy = '1 1 1 1'
+    # bottom is first line
+    subdomain_id = '0 0 1 0
+                    0 1 1 0
+                    0 1 0 0
+                    0 1 1 0'
+  []
+  [add_inner_bdy]
+    type = SideSetsBetweenSubdomainsGenerator
+    input = 'gmg'
+    new_boundary = 'channel'
+    primary_block = '1'
+    paired_block = '0'
+  []
+  [make_nodeset]
+    type = NodeSetsFromSideSetsGenerator
+    input = 'add_inner_bdy'
+    output = true
+  []
+  [follow_channel]
+    type = PolyLineMeshFollowingNodeSetGenerator
+    input = 'make_nodeset'
+    nodeset = 'channel'
+    line_subdomain = '1d_tube'
+
+    # start from the top, going down
+    starting_point = '2 4 0'
+    starting_direction = '0 -1 0'
+    ignore_nodes_behind = true
+
+    # need to hit the channel wall
+    search_radius = '${fparse 1.5}'
+    # Too small and we won't advance, as barycenter won't move
+    # Be careful we should avoid turning back
+    dx = 0.1
+    max_edges = 40
+    num_edges_between_points = 2
+    verbose = true
+  []
+[]

test/tests/meshgenerators/polyline_following_nodeset_mesh_generator/tests

@@ -0,0 +1,15 @@
+[Tests]
+  issues = '#31851'
+  design = 'PolyLineMeshFollowingNodeSetGenerator.md'
+  [2d_channel]
+    type = Exodiff
+    input = '2d_channel.i'
+    cli_args = '--mesh-only'
+    exodiff = '2d_channel_in.e'
+    requirement = 'The system shall be able to generate a polyline following a nodeset.'
+    # mesh generation test
+    recover = false
+    # not implemented for distributed
+    mesh_mode = replicated
+  []
+[]