3D Cutting Capability Implementation

closes idaholab#26579

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

@@ -0,0 +1,68 @@
+# XYZMeshPlaneCutter
+
+!syntax description /Mesh/XYZMeshPlaneCutter
+
+## Overview
+
+The `XYZMeshPlaneCutter` 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/XYZMeshPlaneCutter/input), must be 3D and contain only first-order elements. The cutting plane is specified by [!param](/Mesh/XYZMeshPlaneCutter/plane_normal) and [!param](/Mesh/XYZMeshPlaneCutter/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).
+
+## Methods
+
+`XYZMeshPlaneCutter` first converts all elements of the input mesh into `TET4` elements. Next, the `TET4` elements sliced by the cutting plane are further split into `TET4` elements. 
+
+### Splitting of Non-TET4 Elements
+
+To ensure consistency during cutting, all original mesh elements which are not TET type are first converted into TET elements through splitting. To be specific, each HEX element is split into six TET elements; each PRISM element is split into three TET elements; and each PYRAMID element is split into two TET elements. Details on the splitting approach follow.
+
+!media framework/meshgenerators/hex_split.png
+      style=display: block;margin-left:auto;margin-right:auto;width:32%;float:left;
+      id=hex_split
+      caption=An example of splitting of a HEX element into six TET elements.
+
+!media framework/meshgenerators/prism_split.png
+      style=display: block;margin-left:auto;margin-right:auto;width:32%;float:left;
+      id=prism_split
+      caption=An example of splitting of a PRISM element into three TET elements.
+
+!media framework/meshgenerators/pyramid_split.png
+      style=display: block;margin-left:auto;margin-right:auto;width:36%;float:left;
+      id=pyramid_split
+      caption=An example of splitting of a PYRAMID element into two TET elements.
+
+After this conversion, all the elements become TET elements. In that case, all the subdomain IDs and names can be preserved.
+
+#### Splitting of HEX Elements
+
+There are multiple ways to split one HEX element into multiple TET elements, resulting in either five or six TET elements. Any splitting method must split each of the six quadrilateral faces of a HEX element into two triangles, which can be done in two different ways. As these quadrilateral faces could be shared with neighboring HEX elements, the splitting of the quadrilateral faces on neighboring elements must be performed consistently. To achieve a consistent splitting approach, which will be discussed later in this documentation page, a HEX element needs to be split into six TET elements. An example of this splitting is illustrated in [hex_split]. Note that the splitting approach shown in [hex_split] is not unique and will be discussed later.
+
+#### Splitting of PRISM Elements
+
+A PRISM element can be split into three TET elements. An example of this splitting is illustrated in [prism_split]. Note that the splitting approach shown in [prism_split] is not unique and will be discussed later. Namely, the three quadrilateral faces of a PRISM element need to be split consistently with the neighboring elements.
+
+#### Splitting of PYRAMID Elements
+
+A PYRAMID element can be split into two TET elements. An example of this splitting is illustrated in [pyramid_split]. Note that the splitting approach shown in [pyramid_split] is not unique and will be discussed later. Namely, the one quadrilateral face of a PYRAMID element needs to be split consistently with the neighboring elements.
+
+#### Consistent Splitting for Neighboring Elements
+
+As discussed above, although splitting of non-TET elements into TET elements is not unique, it is crucial to ensure that the splitting of the neighboring elements involves consistent splitting of the quadrilateral faces. To achieve this, the following approach is used. For each quadrilateral face, there are two ways to split it into two triangles, which correspond to the two diagonal lines of the quadrilateral face. Therefore, in order to ensure that one of the two diagonal lines is selected consistently for all the elements, the diagonal line that involves the node with the lowest global node ID among the four nodes of the quadrilateral face is selected.
+
+### Cutting of TET4 Elements along Plane
+
+Once all the elements of the input mesh have been converted into TET elements, the cutting method only needs to be applied to TET elements. First, all the elements that are cut by the given cutting plane are identified. For these involved TET elements, their relationship with the cutting plane can be categorized into one of the six cases shown in [tet_cut]. For each of these six cases, new nodes are created at the intersection points between the cutting plane and the edges of the TET element. Then the red part of the original TET element is removed and new TET element(s) are created as shown in [tet_cut]. The cross-sections created by this cutting procedure are assigned a new boundary ID as defined by [!param](/Mesh/XYZMeshPlaneCutter/cut_face_id)
+
+!media framework/meshgenerators/tet_cut.png
+      style=display: block;margin-left:auto;margin-right:auto;width:75%;
+      id=tet_cut
+      caption=The six possible cases when slicing a TET element. The cutting plane intersection with the element is shown as blue faces. The red part of the original TET element is removed after cutting, while the blue part of the original TET element is kept and split into multiple TET elements if necessary.
+
+## Example Syntax
+
+!listing test/tests/meshgenerators/xyz_mesh_plane_cutter/simple_cut.i block=Mesh/cut
+
+!syntax parameters /Mesh/XYZMeshPlaneCutter
+
+!syntax inputs /Mesh/XYZMeshPlaneCutter
+
+!syntax children /Mesh/XYZMeshPlaneCutter
+

framework/include/meshgenerators/XYZMeshPlaneCutter.h

@@ -0,0 +1,255 @@
+//* 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 "libmesh/face_tri3.h"
+
+/**
+ * This XYZMeshPlaneCutter 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.
+ */
+class XYZMeshPlaneCutter : public MeshGenerator
+{
+public:
+  static InputParameters validParams();
+
+  XYZMeshPlaneCutter(const InputParameters & parameters);
+
+  std::unique_ptr<MeshBase> generate() override;
+
+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.
+  const boundary_id_type _cut_face_id;
+  /// Reference to input mesh pointer
+  std::unique_ptr<MeshBase> & _input;
+
+  /**
+   * Split a HEX8 element into five TET4 element to prepare for the cutting operation.
+   * @param mesh The mesh to be modified
+   * @param elem_id The id of the element to be split
+   * @param converted_elems_ids a vector to record the ids of the newly created TET4 elements
+   */
+  void hexElemSplitter(ReplicatedMesh & mesh,
+                       const dof_id_type elem_id,
+                       std::vector<dof_id_type> & converted_elems_ids);
+
+  /**
+   * Split a PYRAMID5 element into two TET4 element to prepare for the cutting operation.
+   * @param mesh The mesh to be modified
+   * @param elem_id The id of the element to be split
+   * @param converted_elems_ids a vector to record the ids of the newly created TET4 elements
+   */
+  void pyramidElemSplitter(ReplicatedMesh & mesh,
+                           const dof_id_type elem_id,
+                           std::vector<dof_id_type> & converted_elems_ids);
+
+  /**
+   * Split a PRISM6 element into three TET4 element to prepare for the cutting operation.
+   * @param mesh The mesh to be modified
+   * @param elem_id The id of the element to be split
+   * @param converted_elems_ids a vector to record the ids of the newly created TET4 elements
+   */
+  void prismElemSplitter(ReplicatedMesh & mesh,
+                         const dof_id_type elem_id,
+                         std::vector<dof_id_type> & converted_elems_ids);
+
+  /**
+   * 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 index indicating the relation between the point and the plane
+   */
+  short int 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 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 dof_id_type elem_id,
+                      const Point & plane_point,
+                      const Point & plane_normal,
+                      const subdomain_id_type & block_id_to_remove,
+                      std::vector<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
+   */
+  Node * nonDuplicateNodeCreator(ReplicatedMesh & mesh,
+                                 std::vector<Node *> & new_on_plane_nodes,
+                                 const Point new_point);
+
+  /**
+   * Check if two points overlap.
+   * @param point1 The first point
+   * @param point2 The second point
+   * @return true if the two points overlap
+   */
+  bool pointsOverlap(const Point & point1, const Point & point2) const;
+
+  /**
+   * Rotate a HEX8 element nodes to ensure that the node with the minimum id is the first node;
+   * and the node among its three neighboring nodes with the minimum id is the second node.
+   * @param min_id_index The index of the node with the minimum id
+   * @param sec_min_pos The position of the node among its three neighboring nodes with the minimum
+   * id
+   * @param face_rotation A vector to record the rotation of the faces of the HEX8 element
+   */
+  std::vector<unsigned int> nodeRotationHEX8(unsigned int min_id_index,
+                                             unsigned int sec_min_pos,
+                                             std::vector<unsigned int> & face_rotation) const;
+
+  /**
+   * Calculate the three neighboring nodes of a node in a HEX8 element.
+   * @param min_id_index The index of the node with the minimum id
+   * @return a vector of the three neighboring nodes
+   */
+  std::vector<unsigned int> nodeIndicesHEX8(unsigned int min_id_index) const;
+
+  /**
+   * For a rotated nodes that can form a HEX8 element, create a series of four-node set that can
+   * form TET4 elements to replace the original HEX8 element. All the QUAD4 face of the HEX8 element
+   * will be split by the diagonal line that involves the node with the minimum id of that face.
+   * @param hex_nodes A vector of pointers to the nodes that can form a HEX8 element
+   * @param rotated_tet_face_indices A vector of vectors of the original face indices of the HEX8
+   * element corresponding to the faces of the newly created TET4 elements
+   * @return a vector of vectors of pointers to the nodes that can form TET4 elements
+   */
+  std::vector<std::vector<Node *>>
+  hexNodeOptimizer(std::vector<Node *> & hex_nodes,
+                   std::vector<std::vector<unsigned int>> & rotated_tet_face_indices) const;
+
+  /**
+   * For a HEX8 element, determine the direction of the diagonal line of each face that involves the
+   * node with the minimum id of that face.
+   * @param hex_nodes A vector of pointers to the nodes that can form a HEX8 element
+   * @return a vector of boolean values indicating the direction of the diagonal line of each face
+   */
+  std::vector<bool> quadFaceDiagnalDirectionsHex(std::vector<Node *> & hex_nodes) const;
+
+  /**
+   * For a QUAD4 element, determine the direction of the diagonal line that involves the node with
+   * the minimum id of that element.
+   * @param quad_nodes A vector of pointers to the nodes that can form a QUAD4 element
+   * @return a boolean value indicating the direction of the diagonal line
+   */
+  bool quadFaceDiagnalDirection(std::vector<Node *> & quad_nodes) const;
+
+  /**
+   * Creates sets of four nodes indices that can form TET4 elements to replace the original HEX8
+   * element.
+   * @param diagnal_directions A vector of boolean values indicating the direction of the diagonal
+   * line of each face
+   * @param tet_face_indices A vector of vectors of the original face indices of the HEX8 element
+   * corresponding to the faces of the newly created TET4 elements
+   * @return a vector of vectors of node indices that can form TET4 elements
+   */
+  std::vector<std::vector<unsigned int>>
+  tetNodesForHex(const std::vector<bool> diagnal_directions,
+                 std::vector<std::vector<unsigned int>> & tet_face_indices) const;
+
+  /**
+   * Rotate a PRISM6 element nodes to ensure that the node with the minimum id is the first node.
+   * @param min_id_index The index of the node with the minimum id
+   * @param face_rotation A vector to record the rotation of the faces of the PRISM6 element
+   * @return a vector of node indices that can form a PRISM6 element
+   */
+  std::vector<unsigned int> nodeRotationPRISM6(unsigned int min_id_index,
+                                               std::vector<unsigned int> & face_rotation) const;
+
+  /**
+   * For a rotated nodes that can form a PRISM6 element, create a series of four-node set that can
+   * form TET4 elements to replace the original PRISM6 element. All the QUAD4 face of the PRISM6
+   * element will be split by the diagonal line that involves the node with the minimum id of that
+   * face.
+   * @param prism_nodes A vector of pointers to the nodes that can form a PRISM6 element
+   * @param rotated_tet_face_indices A vector of vectors of the original face indices of the PRISM6
+   * element corresponding to the faces of the newly created TET4 elements
+   * @return a vector of vectors of pointers to the nodes that can form TET4 elements
+   */
+  std::vector<std::vector<Node *>>
+  prismNodeOptimizer(std::vector<Node *> & prism_nodes,
+                     std::vector<std::vector<unsigned int>> & rotated_tet_face_indices) const;
+
+  /**
+   * Creates sets of four nodes indices that can form TET4 elements to replace the original PRISM6
+   * element.
+   * @param diagnal_direction A boolean value indicating the direction of the diagonal line of Face
+   * 2
+   * @param tet_face_indices A vector of vectors of the original face indices of the PRISM6 element
+   * corresponding to the faces of the newly created TET4 elements
+   * @return a vector of vectors of node indices that can form TET4 elements
+   */
+  std::vector<std::vector<unsigned int>>
+  tetNodesForPrism(const bool diagnal_direction,
+                   std::vector<std::vector<unsigned int>> & tet_face_indices) const;
+
+  /**
+   * Rotate a PYRAMID5 element nodes to ensure that the node with the minimum id is the first node
+   * for the bottom face.
+   * @param min_id_index The index of the node with the minimum id for the bottom face
+   * @param face_rotation A vector to record the rotation of the faces of the PYRAMID5 element
+   * @return a vector of node indices that can form a PYRAMID5 element
+   */
+  std::vector<unsigned int> nodeRotationPYRAMIND5(unsigned int min_id_index,
+                                                  std::vector<unsigned int> & face_rotation) const;
+
+  /**
+   * For a rotated nodes that can form a PYRAMID5 element, create a series of four-node set that can
+   * form TET4 elements to replace the original PYRAMID5 element. The QUAD4 face of the
+   * PYRAMID5 element will be split by the diagonal line that involves the node with the minimum id
+   * of that face.
+   * @param pyramid_nodes A vector of pointers to the nodes that can form a PYRAMID5 element
+   * @param rotated_tet_face_indices A vector of vectors of the original face indices of the
+   * PYRAMID5 element corresponding to the faces of the newly created TET4 elements
+   * @return a vector of vectors of pointers to the nodes that can form TET4 elements
+   */
+  std::vector<std::vector<Node *>>
+  pyramidNodeOptimizer(std::vector<Node *> & pyramid_nodes,
+                       std::vector<std::vector<unsigned int>> & rotated_tet_face_indices) const;
+};