Add Level Set Based Mesh Cutter

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

@@ -0,0 +1,24 @@
+# CutMeshByLevelSetGenerator
+
+!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)=0`. 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)`). As each element that is across the cutting level set is cut based on the interception points, this mesh generator ensures a smooth cut instead of 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.
+
+## 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/doc/content/source/meshgenerators/CutMeshByPlaneGenerator.md

@@ -4,7 +4,7 @@
 
 ## Overview
 
-The `CutMeshByPlaneGenerator` is basically the 3D version of [`XYMeshLineCutter`](/XYMeshLineCutter.md). It is used to slice a 3D input mesh along a given plane, and discard the portion of the mesh on one side of the plane. The input mesh, given by [!param](/Mesh/CutMeshByPlaneGenerator/input), must be 3D and contain only first-order elements. The cutting plane is specified by [!param](/Mesh/CutMeshByPlaneGenerator/plane_normal) and [!param](/Mesh/CutMeshByPlaneGenerator/plane_point), which are two `libMesh::Point` type data that represent the normal vector of the cutting plane and a point on the cutting plane, respectively. This mesh generator removes the part of the mesh located on the side of the plane in the direction of the normal vector. The mesh is then smoothed to ensure a straight cut instead of a "zigzag" cut along element boundaries as generated by [`PlaneDeletionGenerator`](/PlaneDeletionGenerator.md).
+The `CutMeshByPlaneGenerator` is basically the 3D version of [`XYMeshLineCutter`](/XYMeshLineCutter.md). It is used to slice a 3D input mesh along a given plane, and discard the portion of the mesh on one side of the plane. The input mesh, given by [!param](/Mesh/CutMeshByPlaneGenerator/input), must be 3D and contain only first-order elements. The cutting plane is specified by [!param](/Mesh/CutMeshByPlaneGenerator/plane_normal) and [!param](/Mesh/CutMeshByPlaneGenerator/plane_point), which are two `libMesh::Point` type data that represent the normal vector of the cutting plane and a point on the cutting plane, respectively. This mesh generator removes the part of the mesh located on the side of the plane in the direction of the normal vector. As each element that is across the cutting plane is cut based on the interception points, this mesh generator ensures a straight cut instead of a "zigzag" cut along element boundaries as generated by [`PlaneDeletionGenerator`](/PlaneDeletionGenerator.md).
 
 ## Methods
 

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 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 level set and a line segment defined by two
+   * points separated by the level set.
+   * @param point1 The first point of the line segment
+   * @param point2 The second point of the line segment
+   * @return the intersection point of the level set and the line segment
+   */
+  Point pointPairLevelSetInterception(const Point & point1, const Point & point2);
+
+  /**
+   * For a TET4 elements crossed by the level set, keep the part of the element on the retaining
+   * side of the level set 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) const;
+
+  /**
+   * 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);
+}