Add Level Set Based Mesh Cutter

closes #29703

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

@@ -0,0 +1,24 @@
+# CutMeshByPlaneGenerator
+
+!syntax description /Mesh/CutMeshByLevelSetGenerator
+
+## Overview
+
+The `CutMeshByLevelSetGenerator` is an extended version of [`CutMeshByPlaneGenerator`](/CutMeshByPlaneGenerator.md). It is used to trim a 3D input mesh based on a given level set `f(x,y,z)`. The portion of the input mesh that is within the level set (i.e., `f(x,y,z)<=0`) is retained, while the portion of the input mesh that is outside the level set (i.e., `f(x,y,z)>0`) is discarded. The input mesh, given by [!param](/Mesh/CutMeshByLevelSetGenerator/input), must be 3D and contain only first-order elements. The level set is specified by [!param](/Mesh/CutMeshByLevelSetGenerator/level_set), which can be interpreted by `FParser` as a function of `x`, `y`, and `z` (i.e., `f(x,y,z)`). The mesh is then smoothed to avert a "zigzag" cut along element boundaries.
+
+Using this mesh generator, a 3D structured mesh defined by a bounding box (e.g., generated by [`GeneratedMeshGenerator`](/GeneratedMeshGenerator.md)) can be subtracted into a 3D mesh with its shape define by a given level set. Such a 3D mesh could be further used as an input of [`XYZDelaunayGenerator`](/XYZDelaunayGenerator.md).
+
+## Methods
+
+`CutMeshByLevelSetGenerator` first converts all elements of the input mesh into `TET4` elements. Next, the `TET4` elements sliced by the level set are further split into `TET4` elements. This mesh generator uses exact the same algorithm as its sibling mesh generator, [`CutMeshByPlaneGenerator`](/CutMeshByPlaneGenerator.md). At the first-order element level, cutting by a plane and cutting by a level set are the same.
+
+## Example Syntax
+
+!listing test/tests/meshgenerators/cut_mesh_by_level_set_generator/simple_cut.i block=Mesh/lsc
+
+!syntax parameters /Mesh/CutMeshByLevelSetGenerator
+
+!syntax inputs /Mesh/CutMeshByLevelSetGenerator
+
+!syntax children /Mesh/CutMeshByLevelSetGenerator
+

framework/include/meshgenerators/CutMeshByLevelSetGenerator.h

@@ -0,0 +1,29 @@
+//* 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 "CutMeshByLevelSetGeneratorBase.h"
+
+/**
+ * This CutMeshByLevelSetGenerator object is designed to trim the input mesh by removing all the
+ * elements on outside the give level set with special processing on the elements crossed by the
+ * cutting surface to ensure a smooth cross-section. The output mesh only consists of TET4
+ * elements.
+ */
+class CutMeshByLevelSetGenerator : public CutMeshByLevelSetGeneratorBase
+{
+public:
+  static InputParameters validParams();
+
+  CutMeshByLevelSetGenerator(const InputParameters & parameters);
+
+protected:
+  /// The analytic level set function in the form of a string that can be parsed by FParser
+  const std::string _level_set;
+};

framework/include/meshgenerators/CutMeshByLevelSetGeneratorBase.h

@@ -0,0 +1,100 @@
+//* 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"
+
+#include "FunctionParserUtils.h"
+
+/**
+ * This CutMeshByLevelSetGeneratorBase object is designed to be the base class of mesh generator
+ * that cuts a 3D mesh based on a analytic level set function. The level set function could be
+ * provided explicitly or indirectly.
+ */
+class CutMeshByLevelSetGeneratorBase : public MeshGenerator, public FunctionParserUtils<false>
+{
+public:
+  static InputParameters validParams();
+
+  CutMeshByLevelSetGeneratorBase(const InputParameters & parameters);
+
+  std::unique_ptr<MeshBase> generate() override;
+
+  /// An enum class for style of input polygon size
+  enum class PointLevelSetRelationIndex : short int
+  {
+    level_set_out_side = 1,
+    level_set_in_side = -1,
+    on_level_set = 0
+  };
+
+protected:
+  /// Name of the input mesh
+  const MeshGeneratorName _input_name;
+  /// The boundary id of the surface generated by the cut.
+  boundary_id_type _cut_face_id;
+  /// The boundary name of the surface generated by the cut.
+  const BoundaryName _cut_face_name;
+  /// Reference to input mesh pointer
+  std::unique_ptr<MeshBase> & _input;
+  /// function parser object describing the level set
+  SymFunctionPtr _func_level_set;
+
+  /**
+   * Evaluate the whether a point is on the level set, inside or outside the level set.
+   * @param point The point at which the level set function is to be evaluated
+   * @return the relation of the point to the level set
+   */
+  PointLevelSetRelationIndex pointLevelSetRelation(const Point & point);
+
+  /**
+   * Calculate the intersection point of a plane and a line segment defined by two points separated
+   * by the plane.
+   * @param point1 The first point of the line segment
+   * @param point2 The second point of the line segment
+   * @return the intersection point of the plane and the line segment
+   */
+  Point pointPairLevelSetInterception(const Point & point1, const Point & point2);
+
+  /**
+   * For a TET4 elements crossed by the cutting plane, keep the part of the element on the retaining
+   * side of the plane using a number of TET4 elements.
+   * @param mesh The mesh to be modified
+   * @param bdry_side_list A list that contains the boundary information of the original mesh
+   * @param elem_id The id of the element to be cut
+   * @param block_id_to_remove The subdomain id of the part of the element to be removed
+   * @param new_on_plane_nodes A vector to record the pointers to the newly created nodes on the
+   * cutting plane
+   */
+  void tet4ElemCutter(ReplicatedMesh & mesh,
+                      const std::vector<libMesh::BoundaryInfo::BCTuple> & bdry_side_list,
+                      const dof_id_type elem_id,
+                      const subdomain_id_type & block_id_to_remove,
+                      std::vector<const Node *> & new_on_plane_nodes);
+
+  /**
+   * Check if a position on a plane has already been used as a node in the mesh. If so, return the
+   * pointer to the existing node. If not, create a new node and return the pointer to the new node.
+   * @param mesh The mesh to be modified
+   * @param new_on_plane_nodes A vector to record the pointers to the newly created nodes on the
+   * cutting plane
+   * @param new_point The position of the potential new node
+   * @return a pointer to the existing node or the newly created node
+   */
+  const Node * nonDuplicateNodeCreator(ReplicatedMesh & mesh,
+                                       std::vector<const Node *> & new_on_plane_nodes,
+                                       const Point & new_point);
+
+  /**
+   * Evaluate the level set function at a given point.
+   * @param point The point at which the level set function is to be evaluated
+   * @return the value of the level set function at the given point
+   */
+  Real levelSetEvaluator(const Point & point);
+};

framework/include/meshgenerators/CutMeshByPlaneGenerator.h

@@ -8,99 +8,23 @@
 //* https://www.gnu.org/licenses/lgpl-2.1.html
 
 #pragma once
-#include "MeshGenerator.h"
+#include "CutMeshByLevelSetGeneratorBase.h"
 
 /**
  * This CutMeshByPlaneGenerator object is designed to trim the input mesh by removing all the
  * elements on one side of a given plane with special processing on the elements crossed by the
- * cutting line to ensure a smooth cross-section. The output mesh only consists of TET4 elements.
+ * cutting plane to ensure a smooth cross-section. The output mesh only consists of TET4 elements.
  */
-class CutMeshByPlaneGenerator : public MeshGenerator
+class CutMeshByPlaneGenerator : public CutMeshByLevelSetGeneratorBase
 {
 public:
   static InputParameters validParams();
 
   CutMeshByPlaneGenerator(const InputParameters & parameters);
 
-  std::unique_ptr<MeshBase> generate() override;
-
-  /// An enum class for style of input polygon size
-  enum class PointPlaneRelationIndex : short int
-  {
-    plane_normal_side = 1,
-    opposite_plane_normal_side = -1,
-    on_plane = 0
-  };
-
 protected:
-  /// Name of the input mesh
-  const MeshGeneratorName _input_name;
   /// A point on the cutting plane
   const Point _plane_point;
   /// A normal vector of the cutting plane
   const Point _plane_normal;
-  /// The boundary id of the face generated by the cut.
-  boundary_id_type _cut_face_id;
-  /// The boundary name of the face generated by the cut.
-  const BoundaryName _cut_face_name;
-  /// Reference to input mesh pointer
-  std::unique_ptr<MeshBase> & _input;
-
-  /**
-   * Determine the relation between a point and a plane.
-   * @param point The point of interest to be determined
-   * @param plane_point A point on the plane of interest
-   * @param plane_normal The normal vector of the plane of interest
-   * @return an enum indicating the relation between the point and the plane
-   */
-  PointPlaneRelationIndex pointPlaneRelation(const Point & point,
-                                             const Point & plane_point,
-                                             const Point & plane_normal) const;
-
-  /**
-   * Calculate the intersection point of a plane and a line segment defined by two points separated
-   * by the plane.
-   * @param point1 The first point of the line segment
-   * @param point2 The second point of the line segment
-   * @param plane_point A point on the plane of interest
-   * @param plane_normal The normal vector of the plane of interest
-   * @return the intersection point of the plane and the line segment
-   */
-  Point pointPairPlaneInterception(const Point & point1,
-                                   const Point & point2,
-                                   const Point & plane_point,
-                                   const Point & plane_normal) const;
-
-  /**
-   * For a TET4 elements crossed by the cutting plane, keep the part of the element on the retaining
-   * side of the plane using a number of TET4 elements.
-   * @param mesh The mesh to be modified
-   * @param bdry_side_list A list that contains the boundary information of the original mesh
-   * @param elem_id The id of the element to be cut
-   * @param plane_point A point on the cutting plane
-   * @param plane_normal The normal vector of the cutting plane
-   * @param block_id_to_remove The subdomain id of the part of the element to be removed
-   * @param new_on_plane_nodes A vector to record the pointers to the newly created nodes on the
-   * cutting plane
-   */
-  void tet4ElemCutter(ReplicatedMesh & mesh,
-                      const std::vector<libMesh::BoundaryInfo::BCTuple> & bdry_side_list,
-                      const dof_id_type elem_id,
-                      const Point & plane_point,
-                      const Point & plane_normal,
-                      const subdomain_id_type & block_id_to_remove,
-                      std::vector<const Node *> & new_on_plane_nodes);
-
-  /**
-   * Check if a position on a plane has already been used as a node in the mesh. If so, return the
-   * pointer to the existing node. If not, create a new node and return the pointer to the new node.
-   * @param mesh The mesh to be modified
-   * @param new_on_plane_nodes A vector to record the pointers to the newly created nodes on the
-   * cutting plane
-   * @param new_point The position of the potential new node
-   * @return a pointer to the existing node or the newly created node
-   */
-  const Node * nonDuplicateNodeCreator(ReplicatedMesh & mesh,
-                                       std::vector<const Node *> & new_on_plane_nodes,
-                                       const Point & new_point);
 };

framework/src/meshgenerators/CutMeshByLevelSetGenerator.C

@@ -0,0 +1,59 @@
+//* 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 "CutMeshByLevelSetGenerator.h"
+
+// C++ includes
+#include <cmath>
+
+registerMooseObject("MooseApp", CutMeshByLevelSetGenerator);
+
+InputParameters
+CutMeshByLevelSetGenerator::validParams()
+{
+  InputParameters params = CutMeshByLevelSetGeneratorBase::validParams();
+
+  params.addRequiredParam<std::string>(
+      "level_set", "Level set used to cut the mesh as a function of x, y, and z.");
+  params.addParam<std::vector<std::string>>(
+      "constant_names", {}, "Vector of constants used in the parsed function");
+  params.addParam<std::vector<std::string>>(
+      "constant_expressions",
+      {},
+      "Vector of values for the constants in constant_names (can be an FParser expression)");
+
+  params.addClassDescription(
+      "This CutMeshByLevelSetGenerator object is designed to trim the input mesh by removing all "
+      "the elements on outside the give level set with special processing on the elements crossed "
+      "by the cutting surface to ensure a smooth cross-section. The output mesh only consists of "
+      "TET4 elements.");
+
+  return params;
+}
+
+CutMeshByLevelSetGenerator::CutMeshByLevelSetGenerator(const InputParameters & parameters)
+  : CutMeshByLevelSetGeneratorBase(parameters), _level_set(getParam<std::string>("level_set"))
+{
+  _func_level_set = std::make_shared<SymFunction>();
+  // set FParser internal feature flags
+  setParserFeatureFlags(_func_level_set);
+  if (isParamValid("constant_names") && isParamValid("constant_expressions"))
+    addFParserConstants(_func_level_set,
+                        getParam<std::vector<std::string>>("constant_names"),
+                        getParam<std::vector<std::string>>("constant_expressions"));
+  if (_func_level_set->Parse(_level_set, "x,y,z") >= 0)
+    mooseError("Invalid function f(x,y,z)\n",
+               _func_level_set,
+               "\nin CutMeshByLevelSetGenerator ",
+               name(),
+               ".\n",
+               _func_level_set->ErrorMsg());
+
+  _func_params.resize(3);
+}