diff --git a/framework/doc/content/source/meshgenerators/ManifoldSubdomainGenerator.md b/framework/doc/content/source/meshgenerators/ManifoldSubdomainGenerator.md
new file mode 100644
index 000000000000..b6fde5d68da9
--- /dev/null
+++ b/framework/doc/content/source/meshgenerators/ManifoldSubdomainGenerator.md
@@ -0,0 +1,44 @@
+# ManifoldSubdomainGenerator
+
+!syntax description /Mesh/ManifoldSubdomainGenerator
+
+## Overview
+
+`ManifoldSubdomainGenerator` assigns a subdomain ID based on whether an element's vertex-average
+point lies inside a closed surface mesh (a manifold). The surface must be watertight: every triangle
+must have exactly three neighboring triangles, all elements must be Tri3, and the mesh must be
+consistently oriented. Vertex-average points that lie on the surface within the configured tolerance
+are treated as inside.
+
+The manifold is supplied as a mesh via the MeshGenerator pipeline, not as a raw file path. To use an
+STL file, load it with [FileMeshGenerator](FileMeshGenerator.md) and apply any needed transforms with
+[TransformGenerator](TransformGenerator.md) before passing the result as
+[!param](/Mesh/ManifoldSubdomainGenerator/manifold).
+
+## Example
+
+The example below loads a cube STL as the manifold, translates it into the mesh domain, and assigns
+subdomain `1` to the elements whose vertex-average points lie inside the closed surface.
+
+!listing test/tests/meshgenerators/manifold_subdomain/basic.i block=Mesh
+
+Transforms are composed by chaining mesh generators. The listing above shows a scale step, a rotate
+step, and a translate step in sequence. Users may also limit reassignment to selected existing blocks
+with [!param](/Mesh/ManifoldSubdomainGenerator/restricted_subdomains), or tag elements that lie
+outside the manifold instead by setting
+[!param](/Mesh/ManifoldSubdomainGenerator/location) to `OUTSIDE`.
+
+[!param](/Mesh/ManifoldSubdomainGenerator/surface_tolerance) is an absolute tolerance. Set it
+relative to the manifold length scale and the coordinate noise produced by the mesh export pipeline.
+If coincident vertices differ by more than this tolerance, the manifold validation step can reject an
+otherwise visually closed surface as non-watertight.
+
+This generator uses `Elem::vertex_average()` as its representative point, not the true geometric
+centroid. That keeps the classification inexpensive and matches the sampling convention used by
+similar subdomain tagging generators.
+
+!syntax parameters /Mesh/ManifoldSubdomainGenerator
+
+!syntax inputs /Mesh/ManifoldSubdomainGenerator
+
+!syntax children /Mesh/ManifoldSubdomainGenerator
diff --git a/framework/include/meshgenerators/ManifoldSubdomainGenerator.h b/framework/include/meshgenerators/ManifoldSubdomainGenerator.h
new file mode 100644
index 000000000000..80047fcbc8bf
--- /dev/null
+++ b/framework/include/meshgenerators/ManifoldSubdomainGenerator.h
@@ -0,0 +1,55 @@
+//* 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"
+#include "MooseEnum.h"
+
+class TriangleManifold;
+
+/**
+ * MeshGenerator for defining a subdomain based on whether element vertex averages lie within a
+ * closed manifold.
+ *
+ * The generator intentionally uses Elem::vertex_average() rather than the true geometric centroid
+ * to mirror the inexpensive point-sampling behavior of other subdomain tagging mesh generators in
+ * MOOSE.
+ */
+class ManifoldSubdomainGenerator : public MeshGenerator
+{
+public:
+  /// Declare the input parameters for STL-based subdomain assignment.
+  static InputParameters validParams();
+
+  /// Construct the mesh generator from user input.
+  ManifoldSubdomainGenerator(const InputParameters & parameters);
+
+  /// Apply STL-based subdomain tagging using element vertex averages as the query points.
+  std::unique_ptr<MeshBase> generate() override;
+
+protected:
+  /// Input mesh to modify in place.
+  std::unique_ptr<MeshBase> & _input;
+
+  /// Surface mesh that defines the closed manifold.
+  std::unique_ptr<MeshBase> & _manifold;
+
+  /// Whether to tag the interior or exterior of the STL manifold.
+  const MooseEnum _location;
+
+  /// Target subdomain identifier to assign.
+  const subdomain_id_type _block_id;
+
+  /// Whether retagging is limited to a subset of existing subdomains.
+  const bool _has_restriction;
+
+  /// Absolute tolerance used by the manifold classifier; choose relative to geometry scale/noise.
+  const Real _surface_tolerance;
+};
diff --git a/framework/include/utils/GeometryUtils.h b/framework/include/utils/GeometryUtils.h
index ffbcd4327259..e433cca80dde 100644
--- a/framework/include/utils/GeometryUtils.h
+++ b/framework/include/utils/GeometryUtils.h
@@ -187,4 +187,34 @@ bool arePointsColinear(const Point & p1, const Point & p2, const Point & p3);
  * @return true if the two line segments intersect, false otherwise
  */
 bool segmentsIntersect(const Point & p1, const Point & p2, const Point & p3, const Point & p4);
+
+/**
+ * Compute the squared distance from a point to a 3-D line segment.
+ * @param[in] point point of interest
+ * @param[in] a first segment endpoint
+ * @param[in] b second segment endpoint
+ * @return squared distance from the point to the segment
+ */
+Real pointSegmentDistanceSq(const Point & point, const Point & a, const Point & b);
+
+/**
+ * Compute the squared distance from a point to a 3-D triangle.
+ * @param[in] point point of interest
+ * @param[in] v0 first triangle vertex
+ * @param[in] v1 second triangle vertex
+ * @param[in] v2 third triangle vertex
+ * @return squared distance from the point to the triangle
+ */
+Real
+pointTriangleDistanceSq(const Point & point, const Point & v0, const Point & v1, const Point & v2);
+
+/**
+ * Compute the signed solid angle subtended by one oriented triangle at the query point.
+ * @param[in] point query point
+ * @param[in] v0 first triangle vertex
+ * @param[in] v1 second triangle vertex
+ * @param[in] v2 third triangle vertex
+ * @return signed solid angle in steradians
+ */
+Real solidAngle(const Point & point, const Point & v0, const Point & v1, const Point & v2);
 } // end of namespace geom_utils
diff --git a/framework/include/utils/TriangleManifold.h b/framework/include/utils/TriangleManifold.h
new file mode 100644
index 000000000000..e2ecff9335d7
--- /dev/null
+++ b/framework/include/utils/TriangleManifold.h
@@ -0,0 +1,120 @@
+//* 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 "MooseTypes.h"
+
+#include "libmesh/bounding_box.h"
+#include "libmesh/face_tri3.h"
+#include "libmesh/mesh_tet_interface.h"
+
+#include <iosfwd>
+#include <memory>
+#include <cstdint>
+#include <string>
+#include <unordered_map>
+#include <vector>
+
+class TriangleManifold
+{
+public:
+  TriangleManifold(MeshBase & mesh, const Real surface_tolerance);
+
+  /**
+   * @return True if the point is inside the manifold or lies on the surface.
+   */
+  bool contains(const Point & point) const;
+
+  /**
+   * @return The manifold bounding box.
+   */
+  const libMesh::BoundingBox & boundingBox() const { return _bounding_box; };
+
+  /**
+   * @return The number of triangles in the loaded manifold.
+   */
+  std::size_t numTriangles() const { return _mesh.n_active_elem(); }
+
+private:
+  /**
+   * Result of intersecting the positive x-direction ray with a triangle.
+   *
+   * `Ambiguous` is returned when the hit is too close to an edge, vertex, or the ray origin. In
+   * that case we abandon parity counting and fall back to the more expensive but robust solid-angle
+   * test.
+   */
+  enum class RayIntersection
+  {
+    Miss,
+    Hit,
+    Ambiguous
+  };
+
+  /// Complete post-parse validation and acceleration-structure setup.
+  void finalize();
+
+  /// Build the yz-plane lookup grid used to accelerate +x ray queries.
+  void buildCandidateGrid();
+
+  /// Cheap global bounding-box rejection for containment queries.
+  bool pointInsideBoundingBox(const Point & point) const;
+
+  /// Detect whether a query point lies on or extremely near the manifold surface.
+  bool pointOnSurface(const Point & point) const;
+
+  /// Intersect a positive x-direction ray with a single triangle.
+  RayIntersection rayIntersectsTriangle(const Point & point, const libMesh::Elem & tri) const;
+
+  /// Robust fallback containment query based on accumulated solid angle.
+  bool containsBySolidAngle(const Point & point) const;
+
+  /// Get the subset of triangles whose yz extents may intersect the query ray.
+  std::vector<dof_id_type> rayCandidates(const Point & point) const;
+
+  /// Number of yz-grid bins in the y direction.
+  std::size_t _num_y_cells = 1;
+
+  /// Number of yz-grid bins in the z direction.
+  std::size_t _num_z_cells = 1;
+
+  /// Minimum global y coordinate used to map query points into yz-grid bins.
+  Real _y_min = 0.0;
+
+  /// Minimum global z coordinate used to map query points into yz-grid bins.
+  Real _z_min = 0.0;
+
+  /// Width of one yz-grid cell in the y direction.
+  Real _y_cell_size = 1.0;
+
+  /// Width of one yz-grid cell in the z direction.
+  Real _z_cell_size = 1.0;
+
+  /// Lookup from packed yz-grid cell index to triangles that could intersect the +x query ray.
+  std::unordered_map<dof_id_type, std::vector<dof_id_type>> _ray_grid;
+
+  MeshBase & _mesh;
+
+  /// Absolute tolerance used throughout validation and geometric classification.
+  const Real _surface_tolerance;
+
+  /// Global bounding box of the transformed manifold.
+  const libMesh::BoundingBox _bounding_box;
+
+  // Pre-built point locator for fast close to surface check
+  const std::unique_ptr<libMesh::PointLocatorBase> _point_locator;
+};
+
+class SurfaceChecker : public libMesh::MeshTetInterface
+{
+public:
+  explicit SurfaceChecker(UnstructuredMesh & mesh) : libMesh::MeshTetInterface(mesh) {}
+  void triangulate() override { mooseError("SurfaceChecker is not meant for triangulation."); }
+  std::string improveAndValidate();
+};
diff --git a/framework/src/meshgenerators/ManifoldSubdomainGenerator.C b/framework/src/meshgenerators/ManifoldSubdomainGenerator.C
new file mode 100644
index 000000000000..783752b413cc
--- /dev/null
+++ b/framework/src/meshgenerators/ManifoldSubdomainGenerator.C
@@ -0,0 +1,124 @@
+//* 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 "ManifoldSubdomainGenerator.h"
+
+#include "CastUniquePointer.h"
+#include "MooseMeshUtils.h"
+#include "TriangleManifold.h"
+
+#include "libmesh/elem.h"
+
+#include <set>
+
+registerMooseObject("MooseApp", ManifoldSubdomainGenerator);
+
+InputParameters
+ManifoldSubdomainGenerator::validParams()
+{
+  // This interface intentionally mirrors other point-sampling subdomain tagging generators.
+  MooseEnum location("INSIDE OUTSIDE", "INSIDE");
+
+  InputParameters params = MeshGenerator::validParams();
+  params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify.");
+  params.addRequiredParam<MeshGeneratorName>("manifold", "The mesh defining a closed manifold.");
+  params.addRequiredParam<subdomain_id_type>("block_id", "Subdomain id to assign.");
+  params.addParam<SubdomainName>("block_name", "Optional subdomain name to assign.");
+  params.addParam<MooseEnum>(
+      "location", location, "Control whether the manifold interior or exterior is tagged.");
+  params.addParam<std::vector<SubdomainName>>("restricted_subdomains",
+                                              "Only reset subdomain ID for given subdomains.");
+  params.addRangeCheckedParam<Real>(
+      "surface_tolerance",
+      1e-10,
+      "surface_tolerance>0",
+      "Absolute geometric tolerance used for manifold validation and near-surface classification. "
+      "Choose this relative to the STL length scale and expected coordinate noise.");
+  params.addClassDescription(
+      "Changes the subdomain ID of elements whose vertex-average point lies inside or outside a "
+      "closed manifold defined by surface mesh.");
+  return params;
+}
+
+ManifoldSubdomainGenerator::ManifoldSubdomainGenerator(const InputParameters & parameters)
+  : MeshGenerator(parameters),
+    _input(getMesh("input")),
+    _manifold(getMesh("manifold")),
+    _location(parameters.get<MooseEnum>("location")),
+    _block_id(parameters.get<subdomain_id_type>("block_id")),
+    _has_restriction(isParamValid("restricted_subdomains")),
+    _surface_tolerance(getParam<Real>("surface_tolerance"))
+{
+}
+
+std::unique_ptr<MeshBase>
+ManifoldSubdomainGenerator::generate()
+{
+  std::unique_ptr<MeshBase> mesh = std::move(_input);
+  std::unique_ptr<MeshBase> manifold_mesh = std::move(_manifold);
+
+  // The manifold describes a volume boundary, so the target mesh must also be volumetric.
+  if (mesh->mesh_dimension() != 3)
+    paramError("input", "Only 3D meshes are supported.");
+
+  // Make sure the mesh is prepared
+  if (!manifold_mesh->is_prepared())
+    manifold_mesh->prepare_for_use();
+  // Must be serialized, for now
+  if (!manifold_mesh->is_serial())
+    paramError("manifold", "Manifold mesh must NOT be distributed.");
+  // The manifold must also be 2D
+  if (*(manifold_mesh->elem_dimensions().begin()) != 2 ||
+      *(manifold_mesh->elem_dimensions().rbegin()) != 2)
+    paramError("manifold", "Only 2D meshes are supported.");
+
+  std::set<SubdomainID> restricted_ids;
+  if (_has_restriction)
+  {
+    // Resolve user block names to ids once before the element loop.
+    const auto & names = getParam<std::vector<SubdomainName>>("restricted_subdomains");
+    for (const auto & name : names)
+    {
+      if (!MooseMeshUtils::hasSubdomainName(*mesh, name))
+        paramError("restricted_subdomains", "The block '", name, "' was not found in the mesh");
+
+      restricted_ids.insert(MooseMeshUtils::getSubdomainID(name, *mesh));
+    }
+  }
+
+  if (MooseMeshUtils::hasSubdomainID(*mesh, _block_id) && mesh->processor_id() == 0)
+    mooseInfo("The requested block_id ",
+              _block_id,
+              " already exists on the input mesh. Elements outside the closed manifold that are "
+              "already assigned to this block will remain unchanged.");
+
+  // Build the classifier up front so the per-element loop only performs point-in-manifold queries.
+  TriangleManifold manifold(*manifold_mesh, _surface_tolerance);
+
+  // On distributed meshes, only locally owned active elements are modified on each rank.
+  for (const auto & elem : mesh->active_local_element_ptr_range())
+  {
+    if (_has_restriction && restricted_ids.count(elem->subdomain_id()) == 0)
+      continue;
+
+    // This generator intentionally samples at the vertex average instead of the true geometric
+    // centroid, matching the inexpensive behavior of other subdomain tagging generators.
+    const bool contains = manifold.contains(elem->vertex_average());
+    if ((contains && _location == "INSIDE") || (!contains && _location == "OUTSIDE"))
+      elem->subdomain_id() = _block_id;
+  }
+
+  // Preserve the optional user-facing subdomain name for the assigned id.
+  if (isParamValid("block_name"))
+    mesh->subdomain_name(_block_id) = getParam<SubdomainName>("block_name");
+
+  // Mark derived mesh data as stale because element subdomain assignments have changed.
+  mesh->unset_has_cached_elem_data();
+  return dynamic_pointer_cast<MeshBase>(mesh);
+}
diff --git a/framework/src/utils/GeometryUtils.C b/framework/src/utils/GeometryUtils.C
index 8688e17289d1..04d32816102e 100644
--- a/framework/src/utils/GeometryUtils.C
+++ b/framework/src/utils/GeometryUtils.C
@@ -10,6 +10,10 @@
 #include "GeometryUtils.h"
 #include "MooseUtils.h"
 
+#include <algorithm>
+#include <cmath>
+#include <limits>
+
 namespace geom_utils
 {
 
@@ -342,4 +346,89 @@ segmentsIntersect(const Point & p1, const Point & p2, const Point & p3, const Po
   else
     return false;
 }
+
+Real
+pointSegmentDistanceSq(const Point & point, const Point & a, const Point & b)
+{
+  const Point ab = b - a;
+  const auto length_sq = ab.norm_sq();
+  if (length_sq <= std::numeric_limits<Real>::epsilon())
+    return (point - a).norm_sq();
+
+  const auto t = std::clamp(((point - a) * ab) / length_sq, 0.0, 1.0);
+  const Point projection = a + t * ab;
+  return (point - projection).norm_sq();
+}
+
+Real
+pointTriangleDistanceSq(const Point & point, const Point & v0, const Point & v1, const Point & v2)
+{
+  const Point ab = v1 - v0;
+  const Point ac = v2 - v0;
+  const Point ap = point - v0;
+  const Real d1 = ab * ap;
+  const Real d2 = ac * ap;
+  if (d1 <= 0.0 && d2 <= 0.0)
+    return (point - v0).norm_sq();
+
+  const Point bp = point - v1;
+  const Real d3 = ab * bp;
+  const Real d4 = ac * bp;
+  if (d3 >= 0.0 && d4 <= d3)
+    return (point - v1).norm_sq();
+
+  const Real vc = d1 * d4 - d3 * d2;
+  if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0)
+  {
+    const Real v = d1 / (d1 - d3);
+    const Point projection = v0 + v * ab;
+    return (point - projection).norm_sq();
+  }
+
+  const Point cp = point - v2;
+  const Real d5 = ab * cp;
+  const Real d6 = ac * cp;
+  if (d6 >= 0.0 && d5 <= d6)
+    return (point - v2).norm_sq();
+
+  const Real vb = d5 * d2 - d1 * d6;
+  if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0)
+  {
+    const Real w = d2 / (d2 - d6);
+    const Point projection = v0 + w * ac;
+    return (point - projection).norm_sq();
+  }
+
+  const Real va = d3 * d6 - d5 * d4;
+  if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0)
+  {
+    const Point bc = v2 - v1;
+    const Real w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
+    const Point projection = v1 + w * bc;
+    return (point - projection).norm_sq();
+  }
+
+  const Real denom = 1.0 / (va + vb + vc);
+  const Real v = vb * denom;
+  const Real w = vc * denom;
+  const Point projection = v0 + ab * v + ac * w;
+  return (point - projection).norm_sq();
+}
+
+Real
+solidAngle(const Point & point, const Point & v0, const Point & v1, const Point & v2)
+{
+  const Point a = v0 - point;
+  const Point b = v1 - point;
+  const Point c = v2 - point;
+
+  const Real la = a.norm();
+  const Real lb = b.norm();
+  const Real lc = c.norm();
+
+  const Real numerator = a * (b.cross(c));
+  const Real denominator = la * lb * lc + (a * b) * lc + (b * c) * la + (c * a) * lb;
+
+  return 2.0 * std::atan2(numerator, denominator);
+}
 } // end namespace geom_utils
diff --git a/framework/src/utils/TriangleManifold.C b/framework/src/utils/TriangleManifold.C
new file mode 100644
index 000000000000..bd104c2e8ba5
--- /dev/null
+++ b/framework/src/utils/TriangleManifold.C
@@ -0,0 +1,300 @@
+//* 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 "TriangleManifold.h"
+
+#include "GeometryUtils.h"
+#include "MooseError.h"
+
+#include "libmesh/mesh_tools.h"
+#include "libmesh/unstructured_mesh.h"
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstdint>
+#include <cstring>
+
+namespace TriangleManifoldUtils
+{
+std::uint64_t
+packCell(const std::size_t y_index, const std::size_t z_index)
+{
+  // The acceleration structure is indexed in yz because all rays travel in the +x direction.
+  return (static_cast<std::uint64_t>(y_index) << 32) | static_cast<std::uint64_t>(z_index);
+}
+
+std::size_t
+cellIndex(const Real value, const Real min_value, const Real cell_size, const std::size_t num_cells)
+{
+  // Clamp to the valid range so near-boundary roundoff does not produce an invalid cell id.
+  if (num_cells <= 1)
+    return 0;
+
+  const auto index = static_cast<long long>(std::floor((value - min_value) / cell_size));
+  if (index < 0)
+    return 0;
+  if (index >= static_cast<long long>(num_cells))
+    return num_cells - 1;
+  return static_cast<std::size_t>(index);
+}
+}
+
+TriangleManifold::TriangleManifold(MeshBase & mesh, const Real surface_tolerance)
+  : _mesh(mesh),
+    _surface_tolerance(surface_tolerance),
+    _bounding_box(MeshTools::create_bounding_box(_mesh)),
+    _point_locator(_mesh.sub_point_locator())
+{
+  mooseAssert(_surface_tolerance > 0.0, "surface_tolerance must be strictly positive.");
+
+  mooseAssert(mesh.is_serial(), "Input manifold mesh must be serialized.");
+  mooseAssert(mesh.mesh_dimension() == 2, "Manifold mesh must be a surface.");
+
+  // Finish topology validation and acceleration-structure setup before queries are allowed.
+  finalize();
+
+  _point_locator->set_close_to_point_tol(_surface_tolerance);
+}
+
+bool
+TriangleManifold::contains(const Point & point) const
+{
+  // Most points are rejected here before we perform any per-triangle work.
+  if (!pointInsideBoundingBox(point))
+    return false;
+
+  // By design, near-surface points count as inside for downstream subdomain tagging.
+  if (pointOnSurface(point))
+    return true;
+
+  // Candidate filtering keeps the parity test from touching every triangle in large surfaces.
+  const auto candidates = rayCandidates(point);
+  unsigned int num_hits = 0;
+
+  for (const auto triangle_index : candidates)
+  {
+    const auto tri = _mesh.elem_ptr(triangle_index);
+    // Triangles wholly behind the ray origin cannot contribute to the parity count.
+    if (tri->loose_bounding_box().max()(0) < point(0) - _surface_tolerance)
+      continue;
+
+    switch (rayIntersectsTriangle(point, *tri))
+    {
+      case RayIntersection::Miss:
+        break;
+      case RayIntersection::Hit:
+        // Standard odd/even counting on a closed surface.
+        ++num_hits;
+        break;
+      case RayIntersection::Ambiguous:
+        // Edge and vertex grazing hits are exactly where parity counting becomes brittle.
+        return containsBySolidAngle(point);
+    }
+  }
+
+  return num_hits % 2;
+}
+
+void
+TriangleManifold::finalize()
+{
+  // Validate that the mesh can be used as a manifold
+  // We use the same logic as determining if the mesh can the tetrahedralized.
+  auto umesh = dynamic_cast<UnstructuredMesh *>(&_mesh);
+  SurfaceChecker checker(*umesh);
+  auto msg = checker.improveAndValidate();
+  if (!msg.empty())
+    mooseError(
+        "The inputted surface mesh cannot be treated as manifold for the following reasons:\n",
+        msg);
+
+  buildCandidateGrid();
+}
+
+void
+TriangleManifold::buildCandidateGrid()
+{
+  // Since every containment ray travels in +x, only y and z are needed to choose candidate
+  // triangles for the parity test.
+  const auto extent_y =
+      std::max(boundingBox().max()(1) - boundingBox().min()(1), _surface_tolerance);
+  const auto extent_z =
+      std::max(boundingBox().max()(2) - boundingBox().min()(2), _surface_tolerance);
+
+  // Choose a roughly square yz grid scaled by aspect ratio so candidate lists stay short.
+  const auto target_cells =
+      std::max<dof_id_type>(1, static_cast<dof_id_type>(std::sqrt(numTriangles())));
+  const auto aspect = std::sqrt(extent_y / extent_z);
+
+  _num_y_cells = std::max<dof_id_type>(
+      1, static_cast<dof_id_type>(std::lround(std::sqrt(target_cells) * aspect)));
+  _num_z_cells = std::max<std::size_t>(
+      1, static_cast<std::size_t>(std::ceil(static_cast<Real>(target_cells) / _num_y_cells)));
+
+  _y_min = boundingBox().min()(1);
+  _z_min = boundingBox().min()(2);
+  _y_cell_size = extent_y / _num_y_cells;
+  _z_cell_size = extent_z / _num_z_cells;
+
+  for (const auto elem : _mesh.active_element_ptr_range())
+  {
+    const auto triangle_index = elem->id();
+    const auto bbox = elem->loose_bounding_box();
+
+    const auto y_start =
+        TriangleManifoldUtils::cellIndex(bbox.min()(1), _y_min, _y_cell_size, _num_y_cells);
+    const auto y_stop =
+        TriangleManifoldUtils::cellIndex(bbox.max()(1), _y_min, _y_cell_size, _num_y_cells);
+    const auto z_start =
+        TriangleManifoldUtils::cellIndex(bbox.min()(2), _z_min, _z_cell_size, _num_z_cells);
+    const auto z_stop =
+        TriangleManifoldUtils::cellIndex(bbox.max()(2), _z_min, _z_cell_size, _num_z_cells);
+
+    // Insert the triangle into every grid cell touched by its yz projection.
+    for (const auto iy : make_range(y_start, y_stop + 1))
+      for (const auto iz : make_range(z_start, z_stop + 1))
+        _ray_grid[TriangleManifoldUtils::packCell(iy, iz)].push_back(triangle_index);
+  }
+}
+
+bool
+TriangleManifold::pointInsideBoundingBox(const Point & point) const
+{
+  // Inflate the global box by the tolerance so near-surface points are not rejected too early.
+  return point(0) >= boundingBox().min()(0) - _surface_tolerance &&
+         point(0) <= boundingBox().max()(0) + _surface_tolerance &&
+         point(1) >= boundingBox().min()(1) - _surface_tolerance &&
+         point(1) <= boundingBox().max()(1) + _surface_tolerance &&
+         point(2) >= boundingBox().min()(2) - _surface_tolerance &&
+         point(2) <= boundingBox().max()(2) + _surface_tolerance;
+}
+
+bool
+TriangleManifold::pointOnSurface(const Point & point) const
+{
+  return (*_point_locator)(point) != nullptr;
+}
+
+TriangleManifold::RayIntersection
+TriangleManifold::rayIntersectsTriangle(const Point & point, const libMesh::Elem & tri) const
+{
+  // This is a fixed-direction Moller-Trumbore style ray/triangle intersection test.
+  static const Point direction(1.0, 0.0, 0.0);
+
+  const Point edge1 = tri.node_ref(1) - tri.node_ref(0);
+  const Point edge2 = tri.node_ref(2) - tri.node_ref(0);
+  const Point h = direction.cross(edge2);
+  const Real determinant = edge1 * h;
+  const Real characteristic_length = std::max(edge1.norm(), edge2.norm());
+
+  if (std::abs(determinant) <= _surface_tolerance * characteristic_length)
+    // Nearly parallel triangles are ignored because they do not provide a stable parity event.
+    return RayIntersection::Miss;
+
+  const Real inv_determinant = 1.0 / determinant;
+  const Point s = point - tri.node_ref(0);
+  const Real u = inv_determinant * (s * h);
+  if (u < 0.0 || u > 1.0)
+    return RayIntersection::Miss;
+
+  const Point q = s.cross(edge1);
+  const Real v = inv_determinant * (direction * q);
+  if (v < 0.0 || u + v > 1.0)
+    return RayIntersection::Miss;
+
+  const Real t = inv_determinant * (edge2 * q);
+  if (t < 0.0)
+    // Intersections behind the ray origin do not contribute to +x parity counting.
+    return RayIntersection::Miss;
+  if (t <= _surface_tolerance)
+    // Hits too close to the ray origin are treated as ambiguous boundary situations.
+    return RayIntersection::Ambiguous;
+
+  const Point intersection = point + t * direction;
+  const auto tolerance_sq = _surface_tolerance * _surface_tolerance;
+  if (geom_utils::pointSegmentDistanceSq(intersection, tri.node_ref(0), tri.node_ref(1)) <=
+          tolerance_sq ||
+      geom_utils::pointSegmentDistanceSq(intersection, tri.node_ref(1), tri.node_ref(2)) <=
+          tolerance_sq ||
+      geom_utils::pointSegmentDistanceSq(intersection, tri.node_ref(2), tri.node_ref(0)) <=
+          tolerance_sq)
+    // Edge and vertex hits are where odd/even parity counting is the least reliable.
+    return RayIntersection::Ambiguous;
+
+  return RayIntersection::Hit;
+}
+
+bool
+TriangleManifold::containsBySolidAngle(const Point & point) const
+{
+  // For a closed oriented mesh, the total solid angle is approximately +-4\pi inside and 0
+  // outside. Taking parity over components handles nested shells cleanly.
+  Real total_angle = 0.0;
+  for (const auto elem : _mesh.active_element_ptr_range())
+    total_angle +=
+        geom_utils::solidAngle(point, elem->node_ref(0), elem->node_ref(1), elem->node_ref(2));
+
+  return std::abs(total_angle) > 2.0 * libMesh::pi;
+}
+
+std::vector<dof_id_type>
+TriangleManifold::rayCandidates(const Point & point) const
+{
+  // If the query is outside the manifold's yz extent, the fixed +x ray cannot hit anything.
+  if (point(1) < boundingBox().min()(1) - _surface_tolerance ||
+      point(1) > boundingBox().max()(1) + _surface_tolerance ||
+      point(2) < boundingBox().min()(2) - _surface_tolerance ||
+      point(2) > boundingBox().max()(2) + _surface_tolerance)
+    return {};
+
+  const auto y_index =
+      TriangleManifoldUtils::cellIndex(point(1), _y_min, _y_cell_size, _num_y_cells);
+  const auto z_index =
+      TriangleManifoldUtils::cellIndex(point(2), _z_min, _z_cell_size, _num_z_cells);
+  const auto it = _ray_grid.find(TriangleManifoldUtils::packCell(y_index, z_index));
+  if (it == _ray_grid.end())
+    return {};
+
+  // Return by value to keep the helper simple and avoid exposing internal storage.
+  return it->second;
+}
+
+std::string
+SurfaceChecker::improveAndValidate()
+{
+  using libMesh::MeshTetInterface;
+  auto result = improve_hull_integrity();
+
+  if (result.empty())
+    return "";
+
+  std::ostringstream err_msg;
+  if (result.count(NON_TRI3))
+    err_msg << "- At least one non-Tri3 element was found.\n" << std::endl;
+  if (result.count(MISSING_NEIGHBOR))
+    err_msg << "- At least one triangle without three neighbors was found.\n" << std::endl;
+  if (result.count(EMPTY_MESH))
+    err_msg << "- The surface mesh was empty\n" << std::endl;
+  if (result.count(MISSING_BACKLINK))
+    err_msg << "- At least one triangle neighbor without a return neighbor link was found.\n";
+  if (result.count(BAD_NEIGHBOR_NODES))
+    err_msg << "- At least one triangle neighbor without expected node links was found.\n";
+  if (result.count(NON_ORIENTED))
+    err_msg << "- At least one triangle neighbor with an inconsistent orientation was found.\n";
+  if (result.count(BAD_NEIGHBOR_LINKS))
+    err_msg << "- At least one triangle neighbor with inconsistent node and neighbor links was "
+               "found.\n";
+  if (result.count(DEGENERATE_ELEMENT))
+    err_msg << "- At least one input triangle is degenerate, with near-zero area relative to the "
+               "manifold.\n";
+  if (result.count(DEGENERATE_MESH))
+    err_msg << "- Mesh is degenerate, with zero thickness in at least one direction.\n";
+  return err_msg.str();
+}
diff --git a/test/tests/meshgenerators/manifold_subdomain/basic.i b/test/tests/meshgenerators/manifold_subdomain/basic.i
new file mode 100644
index 000000000000..7c97d03ec45a
--- /dev/null
+++ b/test/tests/meshgenerators/manifold_subdomain/basic.i
@@ -0,0 +1,42 @@
+expected = '56 8'
+
+!include check.i
+
+[Mesh]
+  [gmg]
+    type = GeneratedMeshGenerator
+    dim = 3
+    nx = 4
+    ny = 4
+    nz = 4
+  []
+
+  [stl]
+    type = FileMeshGenerator
+    file = cube_ascii.stl
+  []
+  [scale]
+    type = TransformGenerator
+    input = stl
+    transform = SCALE
+    vector_value = '1 1 1'
+  []
+  [rotate]
+    type = TransformGenerator
+    input = scale
+    transform = ROTATE
+    vector_value = '0 0 0'
+  []
+  [translate]
+    type = TransformGenerator
+    input = rotate
+    transform = TRANSLATE
+    vector_value = '0.5 0.5 0.5'
+  []
+  [apply]
+    type = ManifoldSubdomainGenerator
+    input = gmg
+    manifold = translate
+    block_id = 1
+  []
+[]
diff --git a/test/tests/meshgenerators/manifold_subdomain/check.i b/test/tests/meshgenerators/manifold_subdomain/check.i
new file mode 100644
index 000000000000..813434299edd
--- /dev/null
+++ b/test/tests/meshgenerators/manifold_subdomain/check.i
@@ -0,0 +1,42 @@
+[VectorPostprocessors]
+  [elem_counter]
+    type = ElementCounterWithID
+    id_name = subdomain_id
+    execute_on = 'initial'
+  []
+  [expected_elem]
+    type = ConstantVectorPostprocessor
+    value = '${expected}'
+    execute_on = 'initial'
+  []
+[]
+
+[Postprocessors]
+  [compare]
+    type = VectorPostprocessorComparison
+    comparison_type = EQUALS
+    vectorpostprocessor_a = elem_counter
+    vector_name_a = nelem
+    vectorpostprocessor_b = expected_elem
+    vector_name_b = value
+    execute_on = 'initial'
+  []
+[]
+
+[UserObjects]
+  [terminator]
+    type = Terminator
+    expression = 'compare < 0.5'
+    error_level = ERROR
+    message = 'Manifold did not produce expected number of elements in each subdomain: ${expected}'
+  []
+[]
+
+[Problem]
+  kernel_coverage_check = false
+  solve = false
+[]
+
+[Executioner]
+  type = Steady
+[]
diff --git a/test/tests/meshgenerators/manifold_subdomain/cube_ascii.stl b/test/tests/meshgenerators/manifold_subdomain/cube_ascii.stl
new file mode 100644
index 000000000000..1f97db73005d
--- /dev/null
+++ b/test/tests/meshgenerators/manifold_subdomain/cube_ascii.stl
@@ -0,0 +1,86 @@
+solid cube
+  facet normal 0 0 -1
+    outer loop
+      vertex -0.25 -0.25 -0.25
+      vertex 0.25 -0.25 -0.25
+      vertex 0.25 0.25 -0.25
+    endloop
+  endfacet
+  facet normal 0 0 -1
+    outer loop
+      vertex -0.25 -0.25 -0.25
+      vertex 0.25 0.25 -0.25
+      vertex -0.25 0.25 -0.25
+    endloop
+  endfacet
+  facet normal 0 0 1
+    outer loop
+      vertex -0.25 -0.25 0.25
+      vertex 0.25 0.25 0.25
+      vertex 0.25 -0.25 0.25
+    endloop
+  endfacet
+  facet normal 0 0 1
+    outer loop
+      vertex -0.25 -0.25 0.25
+      vertex -0.25 0.25 0.25
+      vertex 0.25 0.25 0.25
+    endloop
+  endfacet
+  facet normal 0 -1 0
+    outer loop
+      vertex -0.25 -0.25 -0.25
+      vertex 0.25 -0.25 0.25
+      vertex 0.25 -0.25 -0.25
+    endloop
+  endfacet
+  facet normal 0 -1 0
+    outer loop
+      vertex -0.25 -0.25 -0.25
+      vertex -0.25 -0.25 0.25
+      vertex 0.25 -0.25 0.25
+    endloop
+  endfacet
+  facet normal 0 1 0
+    outer loop
+      vertex -0.25 0.25 -0.25
+      vertex 0.25 0.25 -0.25
+      vertex 0.25 0.25 0.25
+    endloop
+  endfacet
+  facet normal 0 1 0
+    outer loop
+      vertex -0.25 0.25 -0.25
+      vertex 0.25 0.25 0.25
+      vertex -0.25 0.25 0.25
+    endloop
+  endfacet
+  facet normal -1 0 0
+    outer loop
+      vertex -0.25 -0.25 -0.25
+      vertex -0.25 0.25 -0.25
+      vertex -0.25 0.25 0.25
+    endloop
+  endfacet
+  facet normal -1 0 0
+    outer loop
+      vertex -0.25 -0.25 -0.25
+      vertex -0.25 0.25 0.25
+      vertex -0.25 -0.25 0.25
+    endloop
+  endfacet
+  facet normal 1 0 0
+    outer loop
+      vertex 0.25 -0.25 -0.25
+      vertex 0.25 0.25 0.25
+      vertex 0.25 0.25 -0.25
+    endloop
+  endfacet
+  facet normal 1 0 0
+    outer loop
+      vertex 0.25 -0.25 -0.25
+      vertex 0.25 -0.25 0.25
+      vertex 0.25 0.25 0.25
+    endloop
+  endfacet
+endsolid cube
diff --git a/test/tests/meshgenerators/manifold_subdomain/distributed.i b/test/tests/meshgenerators/manifold_subdomain/distributed.i
new file mode 100644
index 000000000000..d34243e894ad
--- /dev/null
+++ b/test/tests/meshgenerators/manifold_subdomain/distributed.i
@@ -0,0 +1,22 @@
+[Mesh]
+  [gmg]
+    type = GeneratedMeshGenerator
+    dim = 3
+    nx = 4
+    ny = 4
+    nz = 4
+  []
+  [manifold]
+    type = DistributedRectilinearMeshGenerator
+    dim = 2
+    nx = 10
+    ny = 10
+  []
+
+  [apply]
+    type = ManifoldSubdomainGenerator
+    input = gmg
+    manifold = manifold
+    block_id = 1
+  []
+[]
diff --git a/test/tests/meshgenerators/manifold_subdomain/invalid_input_dimension.i b/test/tests/meshgenerators/manifold_subdomain/invalid_input_dimension.i
new file mode 100644
index 000000000000..2810c36b827d
--- /dev/null
+++ b/test/tests/meshgenerators/manifold_subdomain/invalid_input_dimension.i
@@ -0,0 +1,18 @@
+[Mesh]
+  [gmg]
+    type = GeneratedMeshGenerator
+    dim = 2
+    nx = 4
+    ny = 4
+  []
+  [stl]
+    type = FileMeshGenerator
+    file = cube_ascii.stl
+  []
+  [apply]
+    type = ManifoldSubdomainGenerator
+    input = gmg
+    manifold = stl
+    block_id = 1
+  []
+[]
diff --git a/test/tests/meshgenerators/manifold_subdomain/invalid_manifold_dimension.i b/test/tests/meshgenerators/manifold_subdomain/invalid_manifold_dimension.i
new file mode 100644
index 000000000000..717c76cafab7
--- /dev/null
+++ b/test/tests/meshgenerators/manifold_subdomain/invalid_manifold_dimension.i
@@ -0,0 +1,20 @@
+[Mesh]
+  [gmg]
+    type = GeneratedMeshGenerator
+    dim = 3
+    nx = 4
+    ny = 4
+    nz = 4
+  []
+  [manifold]
+    type = GeneratedMeshGenerator
+    dim = 3
+  []
+
+  [apply]
+    type = ManifoldSubdomainGenerator
+    input = gmg
+    manifold = manifold
+    block_id = 1
+  []
+[]
diff --git a/test/tests/meshgenerators/manifold_subdomain/non_tri.i b/test/tests/meshgenerators/manifold_subdomain/non_tri.i
new file mode 100644
index 000000000000..ff6f850bbc65
--- /dev/null
+++ b/test/tests/meshgenerators/manifold_subdomain/non_tri.i
@@ -0,0 +1,34 @@
+[Mesh]
+  [gmg]
+    type = GeneratedMeshGenerator
+    dim = 3
+    nx = 4
+    ny = 4
+    nz = 4
+  []
+  [manifold_3D]
+    type = GeneratedMeshGenerator
+    dim = 3
+    xmax = 0.5
+    ymax = 0.5
+    zmax = 0.5
+  []
+  [manifold_block]
+    type = LowerDBlockFromSidesetGenerator
+    input = manifold_3D
+    sidesets = 'left right bottom top front back'
+    new_block_id = '99'
+  []
+  [manifold]
+    type = BlockToMeshConverterGenerator
+    input = manifold_block
+    target_blocks = '99'
+  []
+
+  [apply]
+    type = ManifoldSubdomainGenerator
+    input = gmg
+    manifold = manifold
+    block_id = 1
+  []
+[]
diff --git a/test/tests/meshgenerators/manifold_subdomain/open_cube_ascii.stl b/test/tests/meshgenerators/manifold_subdomain/open_cube_ascii.stl
new file mode 100644
index 000000000000..6ee7f1eff3d4
--- /dev/null
+++ b/test/tests/meshgenerators/manifold_subdomain/open_cube_ascii.stl
@@ -0,0 +1,72 @@
+solid open_cube
+  facet normal 0 0 -1
+    outer loop
+      vertex -0.25 -0.25 -0.25
+      vertex 0.25 -0.25 -0.25
+      vertex 0.25 0.25 -0.25
+    endloop
+  endfacet
+  facet normal 0 0 -1
+    outer loop
+      vertex -0.25 -0.25 -0.25
+      vertex 0.25 0.25 -0.25
+      vertex -0.25 0.25 -0.25
+    endloop
+  endfacet
+  facet normal 0 -1 0
+    outer loop
+      vertex -0.25 -0.25 -0.25
+      vertex 0.25 -0.25 0.25
+      vertex 0.25 -0.25 -0.25
+    endloop
+  endfacet
+  facet normal 0 -1 0
+    outer loop
+      vertex -0.25 -0.25 -0.25
+      vertex -0.25 -0.25 0.25
+      vertex 0.25 -0.25 0.25
+    endloop
+  endfacet
+  facet normal 0 1 0
+    outer loop
+      vertex -0.25 0.25 -0.25
+      vertex 0.25 0.25 -0.25
+      vertex 0.25 0.25 0.25
+    endloop
+  endfacet
+  facet normal 0 1 0
+    outer loop
+      vertex -0.25 0.25 -0.25
+      vertex 0.25 0.25 0.25
+      vertex -0.25 0.25 0.25
+    endloop
+  endfacet
+  facet normal -1 0 0
+    outer loop
+      vertex -0.25 -0.25 -0.25
+      vertex -0.25 0.25 -0.25
+      vertex -0.25 0.25 0.25
+    endloop
+  endfacet
+  facet normal -1 0 0
+    outer loop
+      vertex -0.25 -0.25 -0.25
+      vertex -0.25 0.25 0.25
+      vertex -0.25 -0.25 0.25
+    endloop
+  endfacet
+  facet normal 1 0 0
+    outer loop
+      vertex 0.25 -0.25 -0.25
+      vertex 0.25 0.25 0.25
+      vertex 0.25 0.25 -0.25
+    endloop
+  endfacet
+  facet normal 1 0 0
+    outer loop
+      vertex 0.25 -0.25 -0.25
+      vertex 0.25 -0.25 0.25
+      vertex 0.25 0.25 0.25
+    endloop
+  endfacet
+endsolid open_cube
diff --git a/test/tests/meshgenerators/manifold_subdomain/open_surface.i b/test/tests/meshgenerators/manifold_subdomain/open_surface.i
new file mode 100644
index 000000000000..4fdba6ef85bc
--- /dev/null
+++ b/test/tests/meshgenerators/manifold_subdomain/open_surface.i
@@ -0,0 +1,19 @@
+[Mesh]
+  [gmg]
+    type = GeneratedMeshGenerator
+    dim = 3
+    nx = 4
+    ny = 4
+    nz = 4
+  []
+  [stl]
+    type = FileMeshGenerator
+    file = open_cube_ascii.stl
+  []
+  [apply]
+    type = ManifoldSubdomainGenerator
+    input = gmg
+    manifold = stl
+    block_id = 1
+  []
+[]
diff --git a/test/tests/meshgenerators/manifold_subdomain/restricted.i b/test/tests/meshgenerators/manifold_subdomain/restricted.i
new file mode 100644
index 000000000000..3ad27b6b6f68
--- /dev/null
+++ b/test/tests/meshgenerators/manifold_subdomain/restricted.i
@@ -0,0 +1,38 @@
+expected = '32 4 28'
+
+!include check.i
+
+[Mesh]
+  [gmg]
+    type = GeneratedMeshGenerator
+    dim = 3
+    nx = 4
+    ny = 4
+    nz = 4
+  []
+  [seed]
+    type = SubdomainBoundingBoxGenerator
+    input = gmg
+    bottom_left = '0 0 0'
+    top_right = '0.5 1 1'
+    block_id = 2
+  []
+
+  [stl]
+    type = FileMeshGenerator
+    file = cube_ascii.stl
+  []
+  [translate]
+    type = TransformGenerator
+    input = stl
+    transform = TRANSLATE
+    vector_value = '0.5 0.5 0.5'
+  []
+  [apply]
+    type = ManifoldSubdomainGenerator
+    input = seed
+    manifold = translate
+    restricted_subdomains = '2'
+    block_id = 1
+  []
+[]
diff --git a/test/tests/meshgenerators/manifold_subdomain/tests b/test/tests/meshgenerators/manifold_subdomain/tests
new file mode 100644
index 000000000000..a987dc400c3e
--- /dev/null
+++ b/test/tests/meshgenerators/manifold_subdomain/tests
@@ -0,0 +1,68 @@
+[Tests]
+  design = 'meshgenerators/ManifoldSubdomainGenerator.md'
+  issues = '#32763'
+
+  [stl_manifold]
+    [basic_inside]
+      type = RunApp
+      input = basic.i
+      requirement = 'The system shall assign a subdomain to 3D mesh elements whose vertex-average points lie inside a closed STL manifold.'
+    []
+    [outside]
+      type = RunApp
+      input = basic.i
+      cli_args = 'Mesh/apply/location=OUTSIDE expected="8 56"'
+      requirement = 'The system shall assign a subdomain to 3D mesh elements whose vertex-average points lie outside a closed STL manifold.'
+    []
+    [restricted]
+      type = RunApp
+      input = restricted.i
+      requirement = 'The system shall restrict manifold subdomain assignment to selected existing subdomains.'
+    []
+    [rotated_translated]
+      type = RunApp
+      input = basic.i
+      cli_args = 'Mesh/rotate/vector_value="0 0 45" expected="56 8"'
+      requirement = 'The system shall support rigid-body rotation and translation of the STL containment geometry.'
+    []
+    [scaled]
+      type = RunApp
+      input = basic.i
+      cli_args = 'Mesh/scale/vector_value="2 2 2" expected="64"'
+      requirement = 'The system shall support scaling STL geometry before rotation and translation for unit conversion workflows.'
+    []
+  []
+
+  [errors]
+    [invalid_input_dimension]
+      type = RunException
+      input = invalid_input_dimension.i
+      expect_err = 'Only 3D meshes are supported.'
+      requirement = 'The system shall report an error if manifold subdomain assignment is requested on a non-3D mesh.'
+    []
+    [invalid_manifold_dimension]
+      type = RunException
+      input = invalid_manifold_dimension.i
+      expect_err = 'Only 2D meshes are supported.'
+      requirement = 'The system shall report an error if manifold subdomain assignment is requested using non-2D manifold mesh.'
+    []
+    [open_surface]
+      type = RunException
+      input = open_surface.i
+      expect_err = 'At least one triangle without three neighbors was found.'
+      requirement = 'The system shall report an error if the surface is open or non-manifold.'
+    []
+    [non_tri]
+      type = RunException
+      input = non_tri.i
+      expect_err = "At least one non-Tri3 element was found."
+      requirement = 'The system shall report an error if the manifold has non-TRI elements.'
+    []
+    [distributed]
+      type = RunException
+      input = distributed.i
+      expect_err = "Manifold mesh must NOT be distributed."
+      requirement = 'The system shall report an error if the manifold mesh is not serialized.'
+    []
+  []
+[]
diff --git a/unit/include/TriangleManifoldTest.h b/unit/include/TriangleManifoldTest.h
new file mode 100644
index 000000000000..a009f0b382f7
--- /dev/null
+++ b/unit/include/TriangleManifoldTest.h
@@ -0,0 +1,19 @@
+//* 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 "TriangleManifold.h"
+#include "MooseObjectUnitTest.h"
+
+class TriangleManifoldTest : public MooseObjectUnitTest
+{
+public:
+  TriangleManifoldTest() : MooseObjectUnitTest("MooseUnitApp") {}
+};
diff --git a/unit/src/TriangleManifoldTest.C b/unit/src/TriangleManifoldTest.C
new file mode 100644
index 000000000000..85495c8eee22
--- /dev/null
+++ b/unit/src/TriangleManifoldTest.C
@@ -0,0 +1,247 @@
+//* 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 "TriangleManifoldTest.h"
+#include "MooseUnitUtils.h"
+
+#include "libmesh/face_quad4.h"
+#include "libmesh/face_tri3.h"
+#include "libmesh/replicated_mesh.h"
+
+#include "gtest/gtest.h"
+
+namespace
+{
+
+void
+addTri(libMesh::MeshBase & mesh, libMesh::Node * n0, libMesh::Node * n1, libMesh::Node * n2)
+{
+  libMesh::Elem * e = mesh.add_elem(new libMesh::Tri3);
+  e->set_node(0, n0);
+  e->set_node(1, n1);
+  e->set_node(2, n2);
+}
+
+// Watertight unit cube [0,1]^3 as 12 Tri3 elements with outward normals.
+// Nodes: n[0]=(0,0,0) n[1]=(1,0,0) n[2]=(1,1,0) n[3]=(0,1,0)
+//        n[4]=(0,0,1) n[5]=(1,0,1) n[6]=(1,1,1) n[7]=(0,1,1)
+template <typename MeshClass = libMesh::ReplicatedMesh>
+std::unique_ptr<MeshClass>
+makeUnitCubeMesh(const Parallel::Communicator & comm)
+{
+  auto mesh = std::make_unique<MeshClass>(comm);
+  mesh->set_mesh_dimension(2);
+  mesh->set_spatial_dimension(3);
+
+  libMesh::Node * n[8];
+  n[0] = mesh->add_point(Point(0, 0, 0));
+  n[1] = mesh->add_point(Point(1, 0, 0));
+  n[2] = mesh->add_point(Point(1, 1, 0));
+  n[3] = mesh->add_point(Point(0, 1, 0));
+  n[4] = mesh->add_point(Point(0, 0, 1));
+  n[5] = mesh->add_point(Point(1, 0, 1));
+  n[6] = mesh->add_point(Point(1, 1, 1));
+  n[7] = mesh->add_point(Point(0, 1, 1));
+
+  addTri(*mesh, n[0], n[2], n[1]); // bottom -z
+  addTri(*mesh, n[0], n[3], n[2]);
+  addTri(*mesh, n[4], n[5], n[6]); // top +z
+  addTri(*mesh, n[4], n[6], n[7]);
+  addTri(*mesh, n[0], n[1], n[5]); // front -y
+  addTri(*mesh, n[0], n[5], n[4]);
+  addTri(*mesh, n[2], n[3], n[7]); // back +y
+  addTri(*mesh, n[2], n[7], n[6]);
+  addTri(*mesh, n[0], n[4], n[7]); // left -x
+  addTri(*mesh, n[0], n[7], n[3]);
+  addTri(*mesh, n[1], n[2], n[6]); // right +x
+  addTri(*mesh, n[1], n[6], n[5]);
+
+  mesh->prepare_for_use();
+  return mesh;
+}
+
+template std::unique_ptr<libMesh::ReplicatedMesh>
+makeUnitCubeMesh<libMesh::ReplicatedMesh>(const Parallel::Communicator &);
+template std::unique_ptr<libMesh::DistributedMesh>
+makeUnitCubeMesh<libMesh::DistributedMesh>(const Parallel::Communicator &);
+
+// Same cube with the +z face removed (open surface: 10 triangles, boundary edges exposed).
+std::unique_ptr<libMesh::ReplicatedMesh>
+makeOpenCubeMesh(const Parallel::Communicator & comm)
+{
+  auto mesh = std::make_unique<libMesh::ReplicatedMesh>(comm);
+  mesh->set_mesh_dimension(2);
+  mesh->set_spatial_dimension(3);
+
+  libMesh::Node * n[8];
+  n[0] = mesh->add_point(Point(0, 0, 0));
+  n[1] = mesh->add_point(Point(1, 0, 0));
+  n[2] = mesh->add_point(Point(1, 1, 0));
+  n[3] = mesh->add_point(Point(0, 1, 0));
+  n[4] = mesh->add_point(Point(0, 0, 1));
+  n[5] = mesh->add_point(Point(1, 0, 1));
+  n[6] = mesh->add_point(Point(1, 1, 1));
+  n[7] = mesh->add_point(Point(0, 1, 1));
+
+  addTri(*mesh, n[0], n[2], n[1]); // bottom -z
+  addTri(*mesh, n[0], n[3], n[2]);
+  // top face omitted
+  addTri(*mesh, n[0], n[1], n[5]); // front -y
+  addTri(*mesh, n[0], n[5], n[4]);
+  addTri(*mesh, n[2], n[3], n[7]); // back +y
+  addTri(*mesh, n[2], n[7], n[6]);
+  addTri(*mesh, n[0], n[4], n[7]); // left -x
+  addTri(*mesh, n[0], n[7], n[3]);
+  addTri(*mesh, n[1], n[2], n[6]); // right +x
+  addTri(*mesh, n[1], n[6], n[5]);
+
+  mesh->prepare_for_use();
+  return mesh;
+}
+
+// 2D mesh with a single Quad4 element (disallowed by SurfaceChecker).
+std::unique_ptr<libMesh::ReplicatedMesh>
+makeNonTri3Mesh(const Parallel::Communicator & comm)
+{
+  auto mesh = std::make_unique<libMesh::ReplicatedMesh>(comm);
+  mesh->set_mesh_dimension(2);
+  mesh->set_spatial_dimension(3);
+
+  libMesh::Node * n0 = mesh->add_point(Point(0, 0, 0));
+  libMesh::Node * n1 = mesh->add_point(Point(1, 0, 0));
+  libMesh::Node * n2 = mesh->add_point(Point(1, 1, 0));
+  libMesh::Node * n3 = mesh->add_point(Point(0, 1, 0));
+
+  libMesh::Elem * e = mesh->add_elem(new libMesh::Quad4);
+  e->set_node(0, n0);
+  e->set_node(1, n1);
+  e->set_node(2, n2);
+  e->set_node(3, n3);
+
+  mesh->prepare_for_use();
+  return mesh;
+}
+
+// Empty 2D surface mesh (no elements).
+std::unique_ptr<libMesh::ReplicatedMesh>
+makeEmptyMesh(const Parallel::Communicator & comm)
+{
+  auto mesh = std::make_unique<libMesh::ReplicatedMesh>(comm);
+  mesh->set_mesh_dimension(2);
+  mesh->set_spatial_dimension(3);
+  mesh->prepare_for_use();
+  return mesh;
+}
+
+// Single Tri3 with collinear nodes: zero area, no neighbors, degenerate bounding box.
+std::unique_ptr<libMesh::ReplicatedMesh>
+makeDegenerateMesh(const Parallel::Communicator & comm)
+{
+  auto mesh = std::make_unique<libMesh::ReplicatedMesh>(comm);
+  mesh->set_mesh_dimension(2);
+  mesh->set_spatial_dimension(3);
+
+  libMesh::Node * n0 = mesh->add_point(Point(0, 0, 0));
+  libMesh::Node * n1 = mesh->add_point(Point(1, 0, 0));
+  libMesh::Node * n2 = mesh->add_point(Point(2, 0, 0));
+
+  libMesh::Elem * e = mesh->add_elem(new libMesh::Tri3);
+  e->set_node(0, n0);
+  e->set_node(1, n1);
+  e->set_node(2, n2);
+
+  mesh->prepare_for_use();
+  return mesh;
+}
+
+} // namespace
+
+TEST_F(TriangleManifoldTest, containmentBasic)
+{
+  auto mesh = makeUnitCubeMesh(_app->comm());
+  TriangleManifold manifold(*mesh, 1e-10);
+
+  EXPECT_EQ(manifold.numTriangles(), std::size_t{12});
+  EXPECT_TRUE(manifold.contains(Point(0.5, 0.5, 0.5)));   // center
+  EXPECT_FALSE(manifold.contains(Point(1.5, 0.5, 0.5)));  // outside +x
+  EXPECT_FALSE(manifold.contains(Point(-0.5, 0.5, 0.5))); // outside -x
+}
+
+TEST_F(TriangleManifoldTest, surfacePoint)
+{
+  auto mesh = makeUnitCubeMesh(_app->comm());
+  TriangleManifold manifold(*mesh, 1e-6);
+
+  // A point exactly at a mesh vertex is detected as on-surface and counted as inside.
+  EXPECT_TRUE(manifold.contains(Point(0.0, 0.0, 0.0)));
+  // A point on the interior of a face is also on-surface.
+  EXPECT_TRUE(manifold.contains(Point(0.5, 0.5, 0.0)));
+}
+
+TEST_F(TriangleManifoldTest, boundingBox)
+{
+  auto mesh = makeUnitCubeMesh(_app->comm());
+  TriangleManifold manifold(*mesh, 1e-10);
+
+  EXPECT_NEAR(manifold.boundingBox().min()(0), 0.0, 1e-12);
+  EXPECT_NEAR(manifold.boundingBox().min()(1), 0.0, 1e-12);
+  EXPECT_NEAR(manifold.boundingBox().min()(2), 0.0, 1e-12);
+  EXPECT_NEAR(manifold.boundingBox().max()(0), 1.0, 1e-12);
+  EXPECT_NEAR(manifold.boundingBox().max()(1), 1.0, 1e-12);
+  EXPECT_NEAR(manifold.boundingBox().max()(2), 1.0, 1e-12);
+}
+
+TEST_F(TriangleManifoldTest, emptyMesh)
+{
+  auto mesh = makeEmptyMesh(_app->comm());
+#ifdef NDEBUG
+  EXPECT_MOOSEERROR_MSG_CONTAINS(TriangleManifold(*mesh, 1e-10), "surface mesh was empty");
+#else
+  EXPECT_THROW_MSG_CONTAINS(
+      TriangleManifold(*mesh, 1e-10), std::runtime_error, "Manifold mesh must be a surface.");
+#endif
+}
+
+TEST_F(TriangleManifoldTest, openSurface)
+{
+  auto mesh = makeOpenCubeMesh(_app->comm());
+  EXPECT_MOOSEERROR_MSG_CONTAINS(TriangleManifold(*mesh, 1e-10),
+                                 "triangle without three neighbors");
+}
+
+TEST_F(TriangleManifoldTest, nonTri3Element)
+{
+  auto mesh = makeNonTri3Mesh(_app->comm());
+  EXPECT_MOOSEERROR_MSG_CONTAINS(TriangleManifold(*mesh, 1e-10), "non-Tri3 element");
+}
+
+TEST_F(TriangleManifoldTest, degenerateElement)
+{
+  auto mesh = makeDegenerateMesh(_app->comm());
+  EXPECT_MOOSEERROR_MSG_CONTAINS(TriangleManifold(*mesh, 1e-10), "degenerate");
+}
+
+#ifndef NDEBUG
+TEST_F(TriangleManifoldTest, negativeTolerance)
+{
+  auto mesh = makeUnitCubeMesh(_app->comm());
+  EXPECT_THROW_MSG_CONTAINS(TriangleManifold(*mesh, -1e-10),
+                            std::runtime_error,
+                            "surface_tolerance must be strictly positive");
+}
+
+TEST_F(TriangleManifoldTest, distributedMesh)
+{
+  auto mesh = makeUnitCubeMesh<libMesh::DistributedMesh>(_app->comm());
+  mesh->set_distributed();
+  EXPECT_THROW_MSG_CONTAINS(TriangleManifold(*mesh, 1e-10),
+                            std::runtime_error,
+                            "Input manifold mesh must be serialized.");
+}
+#endif