Implementation of STLSubdomainGenerator

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

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;
+};

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

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();
+};

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);
+}