From 2fc0a3e084a5262ea5d051c80ef40e8055ed9f93 Mon Sep 17 00:00:00 2001 From: Yinbin Miao Date: Fri, 6 Dec 2024 10:16:49 -0600 Subject: [PATCH] Add Level Set Based Mesh Cutter closes #29703 --- .../CutMeshByLevelSetGenerator.md | 24 + .../CutMeshByLevelSetGenerator.h | 29 + .../CutMeshByLevelSetGeneratorBase.h | 100 ++++ .../meshgenerators/CutMeshByPlaneGenerator.h | 82 +-- .../CutMeshByLevelSetGenerator.C | 59 ++ .../CutMeshByLevelSetGeneratorBase.C | 509 ++++++++++++++++++ .../meshgenerators/CutMeshByPlaneGenerator.C | 486 +---------------- .../cut_xyzd.i | 43 ++ .../gold/cut_xyzd_out.csv | 3 + .../gold/simple_hex_cut_in.e | Bin 0 -> 25372 bytes .../simple_cut.i | 26 + .../cut_mesh_by_level_set_generator/tests | 25 + 12 files changed, 846 insertions(+), 540 deletions(-) create mode 100644 framework/doc/content/source/meshgenerators/CutMeshByLevelSetGenerator.md create mode 100644 framework/include/meshgenerators/CutMeshByLevelSetGenerator.h create mode 100644 framework/include/meshgenerators/CutMeshByLevelSetGeneratorBase.h create mode 100644 framework/src/meshgenerators/CutMeshByLevelSetGenerator.C create mode 100644 framework/src/meshgenerators/CutMeshByLevelSetGeneratorBase.C create mode 100644 test/tests/meshgenerators/cut_mesh_by_level_set_generator/cut_xyzd.i create mode 100644 test/tests/meshgenerators/cut_mesh_by_level_set_generator/gold/cut_xyzd_out.csv create mode 100644 test/tests/meshgenerators/cut_mesh_by_level_set_generator/gold/simple_hex_cut_in.e create mode 100644 test/tests/meshgenerators/cut_mesh_by_level_set_generator/simple_cut.i create mode 100644 test/tests/meshgenerators/cut_mesh_by_level_set_generator/tests diff --git a/framework/doc/content/source/meshgenerators/CutMeshByLevelSetGenerator.md b/framework/doc/content/source/meshgenerators/CutMeshByLevelSetGenerator.md new file mode 100644 index 000000000000..80fc13e8133a --- /dev/null +++ b/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 + diff --git a/framework/include/meshgenerators/CutMeshByLevelSetGenerator.h b/framework/include/meshgenerators/CutMeshByLevelSetGenerator.h new file mode 100644 index 000000000000..ec948bf405da --- /dev/null +++ b/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; +}; diff --git a/framework/include/meshgenerators/CutMeshByLevelSetGeneratorBase.h b/framework/include/meshgenerators/CutMeshByLevelSetGeneratorBase.h new file mode 100644 index 000000000000..5848cb324e0a --- /dev/null +++ b/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 +{ +public: + static InputParameters validParams(); + + CutMeshByLevelSetGeneratorBase(const InputParameters & parameters); + + std::unique_ptr 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 & _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 & bdry_side_list, + const dof_id_type elem_id, + const subdomain_id_type & block_id_to_remove, + std::vector & 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 & 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); +}; diff --git a/framework/include/meshgenerators/CutMeshByPlaneGenerator.h b/framework/include/meshgenerators/CutMeshByPlaneGenerator.h index b145300c2c23..7ca4dd6ca480 100644 --- a/framework/include/meshgenerators/CutMeshByPlaneGenerator.h +++ b/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 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 & _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 & 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 & 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 & new_on_plane_nodes, - const Point & new_point); }; diff --git a/framework/src/meshgenerators/CutMeshByLevelSetGenerator.C b/framework/src/meshgenerators/CutMeshByLevelSetGenerator.C new file mode 100644 index 000000000000..57593475c723 --- /dev/null +++ b/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 + +registerMooseObject("MooseApp", CutMeshByLevelSetGenerator); + +InputParameters +CutMeshByLevelSetGenerator::validParams() +{ + InputParameters params = CutMeshByLevelSetGeneratorBase::validParams(); + + params.addRequiredParam( + "level_set", "Level set used to cut the mesh as a function of x, y, and z."); + params.addParam>( + "constant_names", {}, "Vector of constants used in the parsed function"); + params.addParam>( + "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("level_set")) +{ + _func_level_set = std::make_shared(); + // set FParser internal feature flags + setParserFeatureFlags(_func_level_set); + if (isParamValid("constant_names") && isParamValid("constant_expressions")) + addFParserConstants(_func_level_set, + getParam>("constant_names"), + getParam>("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); +} diff --git a/framework/src/meshgenerators/CutMeshByLevelSetGeneratorBase.C b/framework/src/meshgenerators/CutMeshByLevelSetGeneratorBase.C new file mode 100644 index 000000000000..afa965e73201 --- /dev/null +++ b/framework/src/meshgenerators/CutMeshByLevelSetGeneratorBase.C @@ -0,0 +1,509 @@ +//* 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 "CutMeshByLevelSetGeneratorBase.h" +#include "MooseMeshElementConversionUtils.h" +#include "MooseMeshUtils.h" + +#include "libmesh/elem.h" +#include "libmesh/boundary_info.h" +#include "libmesh/mesh_base.h" +#include "libmesh/parallel.h" +#include "libmesh/parallel_algebra.h" +#include "libmesh/cell_tet4.h" + +// C++ includes +#include + +InputParameters +CutMeshByLevelSetGeneratorBase::validParams() +{ + InputParameters params = MeshGenerator::validParams(); + params += FunctionParserUtils::validParams(); + + params.addRequiredParam("input", "The input mesh that needs to be trimmed."); + + params.addParam("cut_face_id", + "The boundary id of the face generated by the cut. An " + "id will be automatically assigned if not provided."); + params.addParam("cut_face_name", + "The boundary name of the face generated by the cut."); + + params.addClassDescription( + "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."); + + return params; +} + +CutMeshByLevelSetGeneratorBase::CutMeshByLevelSetGeneratorBase(const InputParameters & parameters) + : MeshGenerator(parameters), + FunctionParserUtils(parameters), + _input_name(getParam("input")), + _cut_face_name(isParamValid("cut_face_name") ? getParam("cut_face_name") + : BoundaryName()), + _input(getMeshByName(_input_name)) +{ + _cut_face_id = isParamValid("cut_face_id") ? getParam("cut_face_id") : -1; +} + +std::unique_ptr +CutMeshByLevelSetGeneratorBase::generate() +{ + auto replicated_mesh_ptr = dynamic_cast(_input.get()); + if (!replicated_mesh_ptr) + paramError("input", "Input is not a replicated mesh, which is required"); + if (*(replicated_mesh_ptr->elem_dimensions().begin()) != 3 || + *(replicated_mesh_ptr->elem_dimensions().rbegin()) != 3) + paramError( + "input", + "Only 3D meshes are supported. Use XYMeshLineCutter for 2D meshes if applicable. Mixed " + "dimensional meshes are not supported at the moment."); + + ReplicatedMesh & mesh = *replicated_mesh_ptr; + + if (!_cut_face_name.empty()) + { + if (MooseMeshUtils::hasBoundaryName(mesh, _cut_face_name)) + { + const boundary_id_type exist_cut_face_id = + MooseMeshUtils::getBoundaryID(_cut_face_name, mesh); + if (_cut_face_id != -1 && _cut_face_id != exist_cut_face_id) + paramError("cut_face_id", + "The cut face boundary name and id are both provided, but they are inconsistent " + "with an existing boundary name which has a different id."); + else + _cut_face_id = exist_cut_face_id; + } + else + { + if (_cut_face_id == -1) + _cut_face_id = MooseMeshUtils::getNextFreeBoundaryID(mesh); + mesh.get_boundary_info().sideset_name(_cut_face_id) = _cut_face_name; + } + } + else if (_cut_face_id == -1) + { + _cut_face_id = MooseMeshUtils::getNextFreeBoundaryID(mesh); + } + + // Subdomain ID for new utility blocks must be new + std::set subdomain_ids_set; + mesh.subdomain_ids(subdomain_ids_set); + const subdomain_id_type max_subdomain_id = *subdomain_ids_set.rbegin(); + const subdomain_id_type block_id_to_remove = max_subdomain_id + 1; + + // For the boolean value in the pair, true means the element is crossed by the cutting plane + // false means the element is on the remaining side + std::vector> cross_and_remained_elems_pre_convert; + + for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end(); + elem_it++) + { + const unsigned int & n_vertices = (*elem_it)->n_vertices(); + unsigned short elem_vertices_counter(0); + for (unsigned int i = 0; i < n_vertices; i++) + { + // We define elem_vertices_counter in this way so that those elements with one face on the + // plane are also processed to have the cut face boundary assigned. + if (pointLevelSetRelation((*(*elem_it)->node_ptr(i))) != + PointLevelSetRelationIndex::level_set_in_side) + elem_vertices_counter++; + } + if (elem_vertices_counter == n_vertices) + (*elem_it)->subdomain_id() = block_id_to_remove; + else + { + // Check if any elements to be processed are not first order + if ((*elem_it)->default_order() != Order::FIRST) + mooseError("Only first order elements are supported for cutting."); + cross_and_remained_elems_pre_convert.push_back( + std::make_pair((*elem_it)->id(), elem_vertices_counter > 0)); + } + } + + std::vector converted_elems_ids_to_cut; + // Then convert these non-TET4 elements into TET4 elements + MooseMeshElementConversionUtils::convert3DMeshToAllTet4(mesh, + cross_and_remained_elems_pre_convert, + converted_elems_ids_to_cut, + block_id_to_remove, + false); + + std::vector new_on_plane_nodes; + // We build the sideset information now, as we only need the information of the elements before + // cutting + BoundaryInfo & boundary_info = mesh.get_boundary_info(); + const auto bdry_side_list = boundary_info.build_side_list(); + // Cut the TET4 Elements + for (const auto & converted_elems_id_to_cut : converted_elems_ids_to_cut) + { + tet4ElemCutter( + mesh, bdry_side_list, converted_elems_id_to_cut, block_id_to_remove, new_on_plane_nodes); + } + + // Delete the block to remove + for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove); + elem_it != mesh.active_subdomain_elements_end(block_id_to_remove); + elem_it++) + mesh.delete_elem(*elem_it); + + mesh.contract(); + mesh.set_isnt_prepared(); + return std::move(_input); +} + +CutMeshByLevelSetGeneratorBase::PointLevelSetRelationIndex +CutMeshByLevelSetGeneratorBase::pointLevelSetRelation(const Point & point) +{ + const Real level_set_value = levelSetEvaluator(point); + if (MooseUtils::absoluteFuzzyLessThan(level_set_value, 0.0)) + return PointLevelSetRelationIndex::level_set_in_side; + else if (MooseUtils::absoluteFuzzyGreaterThan(level_set_value, 0.0)) + return PointLevelSetRelationIndex::level_set_out_side; + else + return PointLevelSetRelationIndex::on_level_set; +} + +Point +CutMeshByLevelSetGeneratorBase::pointPairLevelSetInterception(const Point & point1, + const Point & point2) +{ + Real dist1 = levelSetEvaluator(point1); + Real dist2 = levelSetEvaluator(point2); + + if (MooseUtils::absoluteFuzzyEqual(dist1, 0.0) || MooseUtils::absoluteFuzzyEqual(dist2, 0.0)) + mooseError("At least one of the two points are on the plane."); + if (MooseUtils::absoluteFuzzyGreaterThan(dist1 * dist2, 0.0)) + mooseError("The two points are on the same side of the plane."); + + Point p1(point1); + Point p2(point2); + Real dist = abs(dist1) + abs(dist2); + Point mid_point; + + while (MooseUtils::absoluteFuzzyGreaterThan(dist, 0.0)) + { + mid_point = 0.5 * (p1 + p2); + const Real dist_mid = levelSetEvaluator(mid_point); + // We do not need Fuzzy here as it will be covered by the while loop + if (dist_mid == 0.0) + dist = 0.0; + else if (dist_mid * dist1 < 0.0) + { + p2 = mid_point; + dist2 = levelSetEvaluator(p2); + dist = abs(dist1) + abs(dist2); + } + else + { + p1 = mid_point; + dist1 = levelSetEvaluator(p1); + dist = abs(dist1) + abs(dist2); + } + } + return mid_point; +} + +const Node * +CutMeshByLevelSetGeneratorBase::nonDuplicateNodeCreator( + ReplicatedMesh & mesh, std::vector & new_on_plane_nodes, const Point & new_point) +{ + for (const auto & new_on_plane_node : new_on_plane_nodes) + { + if (MooseUtils::absoluteFuzzyEqual((*new_on_plane_node - new_point).norm(), 0.0)) + return new_on_plane_node; + } + new_on_plane_nodes.push_back(mesh.add_point(new_point)); + return new_on_plane_nodes.back(); +} + +void +CutMeshByLevelSetGeneratorBase::tet4ElemCutter( + ReplicatedMesh & mesh, + const std::vector & bdry_side_list, + const dof_id_type elem_id, + const subdomain_id_type & block_id_to_remove, + std::vector & new_on_plane_nodes) +{ + // Retrieve boundary information for the mesh + BoundaryInfo & boundary_info = mesh.get_boundary_info(); + // Create a list of sidesets involving the element to be split + // It might be complex to assign the boundary id to the new elements + // In TET4, none of the four faces have the same normal vector + // So we are using the normal vector to identify the faces that + // need to be assigned the same boundary id + std::vector elem_side_normal_list; + elem_side_normal_list.resize(4); + for (const auto i : make_range(4)) + { + auto elem_side = mesh.elem_ptr(elem_id)->side_ptr(i); + elem_side_normal_list[i] = (*elem_side->node_ptr(1) - *elem_side->node_ptr(0)) + .cross(*elem_side->node_ptr(2) - *elem_side->node_ptr(1)) + .unit(); + } + + std::vector> elem_side_list; + MooseMeshElementConversionUtils::elementBoundaryInfoCollector( + bdry_side_list, elem_id, 4, elem_side_list); + + std::vector node_plane_relation(4); + std::vector tet4_nodes(4); + std::vector tet4_nodes_on_plane; + std::vector tet4_nodes_outside_plane; + std::vector tet4_nodes_inside_plane; + // Sort tetrahedral nodes based on their positioning wrt the plane + for (unsigned int i = 0; i < 4; i++) + { + tet4_nodes[i] = mesh.elem_ptr(elem_id)->node_ptr(i); + node_plane_relation[i] = pointLevelSetRelation(*tet4_nodes[i]); + if (node_plane_relation[i] == PointLevelSetRelationIndex::on_level_set) + tet4_nodes_on_plane.push_back(tet4_nodes[i]); + else if (node_plane_relation[i] == PointLevelSetRelationIndex::level_set_out_side) + tet4_nodes_outside_plane.push_back(tet4_nodes[i]); + else + tet4_nodes_inside_plane.push_back(tet4_nodes[i]); + } + std::vector elems_tet4; + bool new_elements_created(false); + // No cutting needed if no nodes are outside the plane + if (tet4_nodes_outside_plane.size() == 0) + { + if (tet4_nodes_on_plane.size() == 3) + { + // Record the element for future cross section boundary assignment + elems_tet4.push_back(mesh.elem_ptr(elem_id)); + } + } + // Remove the element if all the nodes are outside the plane + else if (tet4_nodes_inside_plane.size() == 0) + { + mesh.elem_ptr(elem_id)->subdomain_id() = block_id_to_remove; + if (tet4_nodes_on_plane.size() == 3) + { + // I think the neighboring element will be handled, + // So we do not need to assign the cross section boundary here + } + } + // As we have nodes on both sides, six different scenarios are possible + else + { + new_elements_created = true; + if (tet4_nodes_inside_plane.size() == 1 && tet4_nodes_outside_plane.size() == 3) + { + std::vector new_plane_nodes; + // A smaller TET4 element is created, this solution is unique + for (const auto & tet4_node_outside_plane : tet4_nodes_outside_plane) + { + new_plane_nodes.push_back(nonDuplicateNodeCreator( + mesh, + new_on_plane_nodes, + pointPairLevelSetInterception(*tet4_node_outside_plane, *tet4_nodes_inside_plane[0]))); + } + auto new_elem_tet4 = std::make_unique(); + new_elem_tet4->set_node(0) = const_cast(tet4_nodes_inside_plane[0]); + new_elem_tet4->set_node(1) = const_cast(new_plane_nodes[0]); + new_elem_tet4->set_node(2) = const_cast(new_plane_nodes[1]); + new_elem_tet4->set_node(3) = const_cast(new_plane_nodes[2]); + new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id(); + elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4))); + } + else if (tet4_nodes_inside_plane.size() == 2 && tet4_nodes_outside_plane.size() == 2) + { + std::vector new_plane_nodes; + // 3 smaller TET3 elements are created + for (const auto & tet4_node_outside_plane : tet4_nodes_outside_plane) + { + for (const auto & tet4_node_inside_plane : tet4_nodes_inside_plane) + { + new_plane_nodes.push_back(nonDuplicateNodeCreator( + mesh, + new_on_plane_nodes, + pointPairLevelSetInterception(*tet4_node_outside_plane, *tet4_node_inside_plane))); + } + } + std::vector new_elems_nodes = {tet4_nodes_inside_plane[1], + new_plane_nodes[3], + new_plane_nodes[1], + tet4_nodes_inside_plane[0], + new_plane_nodes[2], + new_plane_nodes[0]}; + std::vector> rotated_tet_face_indices; + std::vector> optimized_node_list; + MooseMeshElementConversionUtils::prismNodesToTetNodesDeterminer( + new_elems_nodes, rotated_tet_face_indices, optimized_node_list); + + for (unsigned int i = 0; i < optimized_node_list.size(); i++) + { + auto new_elem_tet4 = std::make_unique(); + new_elem_tet4->set_node(0) = const_cast(optimized_node_list[i][0]); + new_elem_tet4->set_node(1) = const_cast(optimized_node_list[i][1]); + new_elem_tet4->set_node(2) = const_cast(optimized_node_list[i][2]); + new_elem_tet4->set_node(3) = const_cast(optimized_node_list[i][3]); + new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id(); + elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4))); + } + } + else if (tet4_nodes_inside_plane.size() == 3 && tet4_nodes_outside_plane.size() == 1) + { + std::vector new_plane_nodes; + // 3 smaller Tet4 elements are created + for (const auto & tet4_node_inside_plane : tet4_nodes_inside_plane) + { + new_plane_nodes.push_back(nonDuplicateNodeCreator( + mesh, + new_on_plane_nodes, + pointPairLevelSetInterception(*tet4_node_inside_plane, *tet4_nodes_outside_plane[0]))); + } + std::vector new_elems_nodes = {tet4_nodes_inside_plane[0], + tet4_nodes_inside_plane[1], + tet4_nodes_inside_plane[2], + new_plane_nodes[0], + new_plane_nodes[1], + new_plane_nodes[2]}; + std::vector> rotated_tet_face_indices; + std::vector> optimized_node_list; + MooseMeshElementConversionUtils::prismNodesToTetNodesDeterminer( + new_elems_nodes, rotated_tet_face_indices, optimized_node_list); + + for (unsigned int i = 0; i < optimized_node_list.size(); i++) + { + auto new_elem_tet4 = std::make_unique(); + new_elem_tet4->set_node(0) = const_cast(optimized_node_list[i][0]); + new_elem_tet4->set_node(1) = const_cast(optimized_node_list[i][1]); + new_elem_tet4->set_node(2) = const_cast(optimized_node_list[i][2]); + new_elem_tet4->set_node(3) = const_cast(optimized_node_list[i][3]); + new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id(); + elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4))); + } + } + else if (tet4_nodes_inside_plane.size() == 1 && tet4_nodes_outside_plane.size() == 1) + { + auto new_plane_node = nonDuplicateNodeCreator( + mesh, + new_on_plane_nodes, + pointPairLevelSetInterception(*tet4_nodes_inside_plane[0], *tet4_nodes_outside_plane[0])); + // A smaller Tet4 is created, this solution is unique + auto new_elem_tet4 = std::make_unique(); + new_elem_tet4->set_node(0) = const_cast(new_plane_node); + new_elem_tet4->set_node(1) = const_cast(tet4_nodes_on_plane[0]); + new_elem_tet4->set_node(2) = const_cast(tet4_nodes_on_plane[1]); + new_elem_tet4->set_node(3) = const_cast(tet4_nodes_inside_plane[0]); + new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id(); + elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4))); + } + else if (tet4_nodes_inside_plane.size() == 1 && tet4_nodes_outside_plane.size() == 2) + { + std::vector new_plane_nodes; + // A smaller Tet4 element is created, this solution is unique + for (const auto & tet4_node_outside_plane : tet4_nodes_outside_plane) + { + new_plane_nodes.push_back(nonDuplicateNodeCreator( + mesh, + new_on_plane_nodes, + pointPairLevelSetInterception(*tet4_node_outside_plane, *tet4_nodes_inside_plane[0]))); + } + auto new_elem_tet4 = std::make_unique(); + new_elem_tet4->set_node(0) = const_cast(new_plane_nodes[0]); + new_elem_tet4->set_node(1) = const_cast(new_plane_nodes[1]); + new_elem_tet4->set_node(2) = const_cast(tet4_nodes_on_plane[0]); + new_elem_tet4->set_node(3) = const_cast(tet4_nodes_inside_plane[0]); + new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id(); + elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4))); + } + else if (tet4_nodes_inside_plane.size() == 2 && tet4_nodes_outside_plane.size() == 1) + { + std::vector new_plane_nodes; + // 2 smaller TET4 elements are created + for (const auto & tet4_node_inside_plane : tet4_nodes_inside_plane) + { + new_plane_nodes.push_back(nonDuplicateNodeCreator( + mesh, + new_on_plane_nodes, + pointPairLevelSetInterception(*tet4_node_inside_plane, *tet4_nodes_outside_plane[0]))); + } + std::vector new_elems_nodes = {tet4_nodes_inside_plane[0], + tet4_nodes_inside_plane[1], + new_plane_nodes[1], + new_plane_nodes[0], + tet4_nodes_on_plane[0]}; + std::vector> rotated_tet_face_indices; + std::vector> optimized_node_list; + MooseMeshElementConversionUtils::pyramidNodesToTetNodesDeterminer( + new_elems_nodes, rotated_tet_face_indices, optimized_node_list); + + for (unsigned int i = 0; i < optimized_node_list.size(); i++) + { + auto new_elem_tet4 = std::make_unique(); + new_elem_tet4->set_node(0) = const_cast(optimized_node_list[i][0]); + new_elem_tet4->set_node(1) = const_cast(optimized_node_list[i][1]); + new_elem_tet4->set_node(2) = const_cast(optimized_node_list[i][2]); + new_elem_tet4->set_node(3) = const_cast(optimized_node_list[i][3]); + new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id(); + elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4))); + } + } + else + mooseError("Unexpected scenario."); + + mesh.elem_ptr(elem_id)->subdomain_id() = block_id_to_remove; + } + + for (auto & elem_tet4 : elems_tet4) + { + if (new_elements_created) + { + if (elem_tet4->volume() < 0.0) + { + Node * temp = elem_tet4->node_ptr(0); + elem_tet4->set_node(0) = elem_tet4->node_ptr(1); + elem_tet4->set_node(1) = temp; + } + } + // Find the boundary id of the new element + for (unsigned int i = 0; i < 4; i++) + { + const Point & side_pt_0 = *elem_tet4->side_ptr(i)->node_ptr(0); + const Point & side_pt_1 = *elem_tet4->side_ptr(i)->node_ptr(1); + const Point & side_pt_2 = *elem_tet4->side_ptr(i)->node_ptr(2); + + const Point side_normal = (side_pt_1 - side_pt_0).cross(side_pt_2 - side_pt_1).unit(); + for (unsigned int j = 0; j < 4; j++) + { + if (new_elements_created) + { + if (MooseUtils::absoluteFuzzyEqual(side_normal * elem_side_normal_list[j], 1.0)) + { + for (const auto & side_info : elem_side_list[j]) + { + boundary_info.add_side(elem_tet4, i, side_info); + } + } + } + } + if (MooseUtils::absoluteFuzzyEqual(levelSetEvaluator(side_pt_0), 0.0) && + MooseUtils::absoluteFuzzyEqual(levelSetEvaluator(side_pt_1), 0.0) && + MooseUtils::absoluteFuzzyEqual(levelSetEvaluator(side_pt_2), 0.0)) + { + + boundary_info.add_side(elem_tet4, i, _cut_face_id); + } + } + } +} + +Real +CutMeshByLevelSetGeneratorBase::levelSetEvaluator(const Point & point) +{ + _func_params[0] = point(0); + _func_params[1] = point(1); + _func_params[2] = point(2); + return evaluate(_func_level_set); +} diff --git a/framework/src/meshgenerators/CutMeshByPlaneGenerator.C b/framework/src/meshgenerators/CutMeshByPlaneGenerator.C index 3c82f4f2c5c1..282e5253d1c3 100644 --- a/framework/src/meshgenerators/CutMeshByPlaneGenerator.C +++ b/framework/src/meshgenerators/CutMeshByPlaneGenerator.C @@ -8,15 +8,6 @@ //* https://www.gnu.org/licenses/lgpl-2.1.html #include "CutMeshByPlaneGenerator.h" -#include "MooseMeshElementConversionUtils.h" -#include "MooseMeshUtils.h" - -#include "libmesh/elem.h" -#include "libmesh/boundary_info.h" -#include "libmesh/mesh_base.h" -#include "libmesh/parallel.h" -#include "libmesh/parallel_algebra.h" -#include "libmesh/cell_tet4.h" // C++ includes #include @@ -26,471 +17,44 @@ registerMooseObject("MooseApp", CutMeshByPlaneGenerator); InputParameters CutMeshByPlaneGenerator::validParams() { - InputParameters params = MeshGenerator::validParams(); - - params.addRequiredParam("input", "The input mesh that needs to be trimmed."); + InputParameters params = CutMeshByLevelSetGeneratorBase::validParams(); params.addRequiredParam("plane_point", "A point on the plane."); params.addRequiredParam("plane_normal", "The normal vector of the plane."); - params.addParam("cut_face_id", - "The boundary id of the face generated by the cut. An " - "id will be automatically assigned if not provided."); - params.addParam("cut_face_name", - "The boundary name of the face generated by the cut."); params.addClassDescription( "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 " + "the cutting plane to ensure a smooth cross-section. The output mesh only consists of TET4 " "elements."); return params; } CutMeshByPlaneGenerator::CutMeshByPlaneGenerator(const InputParameters & parameters) - : MeshGenerator(parameters), - _input_name(getParam("input")), + : CutMeshByLevelSetGeneratorBase(parameters), _plane_point(getParam("plane_point")), - _plane_normal(getParam("plane_normal").unit()), - _cut_face_name(isParamValid("cut_face_name") ? getParam("cut_face_name") - : BoundaryName()), - _input(getMeshByName(_input_name)) -{ - _cut_face_id = isParamValid("cut_face_id") ? getParam("cut_face_id") : -1; -} - -std::unique_ptr -CutMeshByPlaneGenerator::generate() -{ - auto replicated_mesh_ptr = dynamic_cast(_input.get()); - if (!replicated_mesh_ptr) - paramError("input", "Input is not a replicated mesh, which is required"); - if (*(replicated_mesh_ptr->elem_dimensions().begin()) != 3 || - *(replicated_mesh_ptr->elem_dimensions().rbegin()) != 3) - paramError("input", - "Only 3D meshes are supported. Use XYMeshLineCutter for 2D meshes. Mixed " - "dimensional meshes are not supported at the moment."); - - ReplicatedMesh & mesh = *replicated_mesh_ptr; - - if (!_cut_face_name.empty()) - { - if (MooseMeshUtils::hasBoundaryName(mesh, _cut_face_name)) - { - const boundary_id_type exist_cut_face_id = - MooseMeshUtils::getBoundaryID(_cut_face_name, mesh); - if (_cut_face_id != -1 && _cut_face_id != exist_cut_face_id) - paramError("cut_face_id", - "The cut face boundary name and id are both provided, but they are inconsistent " - "with an existing boundary name which has a different id."); - else - _cut_face_id = exist_cut_face_id; - } - else - { - if (_cut_face_id == -1) - _cut_face_id = MooseMeshUtils::getNextFreeBoundaryID(mesh); - mesh.get_boundary_info().sideset_name(_cut_face_id) = _cut_face_name; - } - } - else if (_cut_face_id == -1) - { - _cut_face_id = MooseMeshUtils::getNextFreeBoundaryID(mesh); - } - - // Subdomain ID for new utility blocks must be new - std::set subdomain_ids_set; - mesh.subdomain_ids(subdomain_ids_set); - const subdomain_id_type max_subdomain_id = *subdomain_ids_set.rbegin(); - const subdomain_id_type block_id_to_remove = max_subdomain_id + 1; - - // For the boolean value in the pair, true means the element is crossed by the cutting plane - // false means the element is on the remaining side - std::vector> cross_and_remained_elems_pre_convert; - - for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end(); - elem_it++) - { - const unsigned int & n_vertices = (*elem_it)->n_vertices(); - unsigned short elem_vertices_counter(0); - for (unsigned int i = 0; i < n_vertices; i++) - { - // We define elem_vertices_counter in this way so that those elements with one face on the - // plane are also processed to have the cut face boundary assigned. - if (pointPlaneRelation((*(*elem_it)->node_ptr(i)), _plane_point, _plane_normal) != - PointPlaneRelationIndex::opposite_plane_normal_side) - elem_vertices_counter++; - } - if (elem_vertices_counter == n_vertices) - (*elem_it)->subdomain_id() = block_id_to_remove; - else - { - // Check if any elements to be processed are not first order - if ((*elem_it)->default_order() != Order::FIRST) - mooseError("Only first order elements are supported for cutting."); - cross_and_remained_elems_pre_convert.push_back( - std::make_pair((*elem_it)->id(), elem_vertices_counter > 0)); - } - } - - std::vector converted_elems_ids_to_cut; - // Then convert these non-TET4 elements into TET4 elements - MooseMeshElementConversionUtils::convert3DMeshToAllTet4(mesh, - cross_and_remained_elems_pre_convert, - converted_elems_ids_to_cut, - block_id_to_remove, - false); - - std::vector new_on_plane_nodes; - // We build the sideset information now, as we only need the information of the elements before - // cutting - BoundaryInfo & boundary_info = mesh.get_boundary_info(); - const auto bdry_side_list = boundary_info.build_side_list(); - // Cut the TET4 Elements - for (const auto & converted_elems_id_to_cut : converted_elems_ids_to_cut) - { - tet4ElemCutter(mesh, - bdry_side_list, - converted_elems_id_to_cut, - _plane_point, - _plane_normal, - block_id_to_remove, - new_on_plane_nodes); - } - - // Delete the block to remove - for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove); - elem_it != mesh.active_subdomain_elements_end(block_id_to_remove); - elem_it++) - mesh.delete_elem(*elem_it); - - mesh.contract(); - mesh.set_isnt_prepared(); - return std::move(_input); -} - -CutMeshByPlaneGenerator::PointPlaneRelationIndex -CutMeshByPlaneGenerator::pointPlaneRelation(const Point & point, - const Point & plane_point, - const Point & plane_normal) const + _plane_normal(getParam("plane_normal").unit()) { - const Real dist = plane_normal * (point - plane_point); - if (MooseUtils::absoluteFuzzyLessThan(dist, 0.0)) - return PointPlaneRelationIndex::opposite_plane_normal_side; - else if (MooseUtils::absoluteFuzzyGreaterThan(dist, 0.0)) - return PointPlaneRelationIndex::plane_normal_side; - else - return PointPlaneRelationIndex::on_plane; -} - -Point -CutMeshByPlaneGenerator::pointPairPlaneInterception(const Point & point1, - const Point & point2, - const Point & plane_point, - const Point & plane_normal) const -{ - const Real dist1 = plane_normal * (point1 - plane_point); - const Real dist2 = plane_normal * (point2 - plane_point); - - if (MooseUtils::absoluteFuzzyEqual(dist1, 0.0) || MooseUtils::absoluteFuzzyEqual(dist2, 0.0)) - mooseError("At least one of the two points are on the plane."); - if (MooseUtils::absoluteFuzzyGreaterThan(dist1 * dist2, 0.0)) - mooseError("The two points are on the same side of the plane."); - - const Real dist_ratio = dist1 / (dist1 - dist2); - return point1 + dist_ratio * (point2 - point1); -} - -const Node * -CutMeshByPlaneGenerator::nonDuplicateNodeCreator(ReplicatedMesh & mesh, - std::vector & new_on_plane_nodes, - const Point & new_point) -{ - for (const auto & new_on_plane_node : new_on_plane_nodes) - { - if (MooseUtils::absoluteFuzzyEqual((*new_on_plane_node - new_point).norm(), 0.0)) - return new_on_plane_node; - } - new_on_plane_nodes.push_back(mesh.add_point(new_point)); - return new_on_plane_nodes.back(); -} - -void -CutMeshByPlaneGenerator::tet4ElemCutter( - ReplicatedMesh & mesh, - const std::vector & 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 & new_on_plane_nodes) -{ - // Retrieve boundary information for the mesh - BoundaryInfo & boundary_info = mesh.get_boundary_info(); - // Create a list of sidesets involving the element to be split - // It might be complex to assign the boundary id to the new elements - // In TET4, none of the four faces have the same normal vector - // So we are using the normal vector to identify the faces that - // need to be assigned the same boundary id - std::vector elem_side_normal_list; - elem_side_normal_list.resize(4); - for (const auto i : make_range(4)) - { - auto elem_side = mesh.elem_ptr(elem_id)->side_ptr(i); - elem_side_normal_list[i] = (*elem_side->node_ptr(1) - *elem_side->node_ptr(0)) - .cross(*elem_side->node_ptr(2) - *elem_side->node_ptr(1)) - .unit(); - } - - std::vector> elem_side_list; - MooseMeshElementConversionUtils::elementBoundaryInfoCollector( - bdry_side_list, elem_id, 4, elem_side_list); - - std::vector node_plane_relation(4); - std::vector tet4_nodes(4); - std::vector tet4_nodes_on_plane; - std::vector tet4_nodes_outside_plane; - std::vector tet4_nodes_inside_plane; - // Sort tetrahedral nodes based on their positioning wrt the plane - for (unsigned int i = 0; i < 4; i++) - { - tet4_nodes[i] = mesh.elem_ptr(elem_id)->node_ptr(i); - node_plane_relation[i] = pointPlaneRelation(*tet4_nodes[i], plane_point, plane_normal); - if (node_plane_relation[i] == PointPlaneRelationIndex::on_plane) - tet4_nodes_on_plane.push_back(tet4_nodes[i]); - else if (node_plane_relation[i] == PointPlaneRelationIndex::plane_normal_side) - tet4_nodes_outside_plane.push_back(tet4_nodes[i]); - else - tet4_nodes_inside_plane.push_back(tet4_nodes[i]); - } - std::vector elems_tet4; - bool new_elements_created(false); - // No cutting needed if no nodes are outside the plane - if (tet4_nodes_outside_plane.size() == 0) - { - if (tet4_nodes_on_plane.size() == 3) - { - // Record the element for future cross section boundary assignment - elems_tet4.push_back(mesh.elem_ptr(elem_id)); - } - } - // Remove the element if all the nodes are outside the plane - else if (tet4_nodes_inside_plane.size() == 0) - { - mesh.elem_ptr(elem_id)->subdomain_id() = block_id_to_remove; - if (tet4_nodes_on_plane.size() == 3) - { - // I think the neighboring element will be handled, - // So we do not need to assign the cross section boundary here - } - } - // As we have nodes on both sides, six different scenarios are possible - else - { - new_elements_created = true; - if (tet4_nodes_inside_plane.size() == 1 && tet4_nodes_outside_plane.size() == 3) - { - std::vector new_plane_nodes; - // A smaller TET4 element is created, this solution is unique - for (const auto & tet4_node_outside_plane : tet4_nodes_outside_plane) - { - new_plane_nodes.push_back(nonDuplicateNodeCreator( - mesh, - new_on_plane_nodes, - pointPairPlaneInterception( - *tet4_node_outside_plane, *tet4_nodes_inside_plane[0], plane_point, plane_normal))); - } - auto new_elem_tet4 = std::make_unique(); - new_elem_tet4->set_node(0) = const_cast(tet4_nodes_inside_plane[0]); - new_elem_tet4->set_node(1) = const_cast(new_plane_nodes[0]); - new_elem_tet4->set_node(2) = const_cast(new_plane_nodes[1]); - new_elem_tet4->set_node(3) = const_cast(new_plane_nodes[2]); - new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id(); - elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4))); - } - else if (tet4_nodes_inside_plane.size() == 2 && tet4_nodes_outside_plane.size() == 2) - { - std::vector new_plane_nodes; - // 3 smaller TET3 elements are created - for (const auto & tet4_node_outside_plane : tet4_nodes_outside_plane) - { - for (const auto & tet4_node_inside_plane : tet4_nodes_inside_plane) - { - new_plane_nodes.push_back(nonDuplicateNodeCreator( - mesh, - new_on_plane_nodes, - pointPairPlaneInterception( - *tet4_node_outside_plane, *tet4_node_inside_plane, plane_point, plane_normal))); - } - } - std::vector new_elems_nodes = {tet4_nodes_inside_plane[1], - new_plane_nodes[3], - new_plane_nodes[1], - tet4_nodes_inside_plane[0], - new_plane_nodes[2], - new_plane_nodes[0]}; - std::vector> rotated_tet_face_indices; - std::vector> optimized_node_list; - MooseMeshElementConversionUtils::prismNodesToTetNodesDeterminer( - new_elems_nodes, rotated_tet_face_indices, optimized_node_list); - - for (unsigned int i = 0; i < optimized_node_list.size(); i++) - { - auto new_elem_tet4 = std::make_unique(); - new_elem_tet4->set_node(0) = const_cast(optimized_node_list[i][0]); - new_elem_tet4->set_node(1) = const_cast(optimized_node_list[i][1]); - new_elem_tet4->set_node(2) = const_cast(optimized_node_list[i][2]); - new_elem_tet4->set_node(3) = const_cast(optimized_node_list[i][3]); - new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id(); - elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4))); - } - } - else if (tet4_nodes_inside_plane.size() == 3 && tet4_nodes_outside_plane.size() == 1) - { - std::vector new_plane_nodes; - // 3 smaller Tet4 elements are created - for (const auto & tet4_node_inside_plane : tet4_nodes_inside_plane) - { - new_plane_nodes.push_back(nonDuplicateNodeCreator( - mesh, - new_on_plane_nodes, - pointPairPlaneInterception( - *tet4_node_inside_plane, *tet4_nodes_outside_plane[0], plane_point, plane_normal))); - } - std::vector new_elems_nodes = {tet4_nodes_inside_plane[0], - tet4_nodes_inside_plane[1], - tet4_nodes_inside_plane[2], - new_plane_nodes[0], - new_plane_nodes[1], - new_plane_nodes[2]}; - std::vector> rotated_tet_face_indices; - std::vector> optimized_node_list; - MooseMeshElementConversionUtils::prismNodesToTetNodesDeterminer( - new_elems_nodes, rotated_tet_face_indices, optimized_node_list); - - for (unsigned int i = 0; i < optimized_node_list.size(); i++) - { - auto new_elem_tet4 = std::make_unique(); - new_elem_tet4->set_node(0) = const_cast(optimized_node_list[i][0]); - new_elem_tet4->set_node(1) = const_cast(optimized_node_list[i][1]); - new_elem_tet4->set_node(2) = const_cast(optimized_node_list[i][2]); - new_elem_tet4->set_node(3) = const_cast(optimized_node_list[i][3]); - new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id(); - elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4))); - } - } - else if (tet4_nodes_inside_plane.size() == 1 && tet4_nodes_outside_plane.size() == 1) - { - auto new_plane_node = - nonDuplicateNodeCreator(mesh, - new_on_plane_nodes, - pointPairPlaneInterception(*tet4_nodes_inside_plane[0], - *tet4_nodes_outside_plane[0], - plane_point, - plane_normal)); - // A smaller Tet4 is created, this solution is unique - auto new_elem_tet4 = std::make_unique(); - new_elem_tet4->set_node(0) = const_cast(new_plane_node); - new_elem_tet4->set_node(1) = const_cast(tet4_nodes_on_plane[0]); - new_elem_tet4->set_node(2) = const_cast(tet4_nodes_on_plane[1]); - new_elem_tet4->set_node(3) = const_cast(tet4_nodes_inside_plane[0]); - new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id(); - elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4))); - } - else if (tet4_nodes_inside_plane.size() == 1 && tet4_nodes_outside_plane.size() == 2) - { - std::vector new_plane_nodes; - // A smaller Tet4 element is created, this solution is unique - for (const auto & tet4_node_outside_plane : tet4_nodes_outside_plane) - { - new_plane_nodes.push_back(nonDuplicateNodeCreator( - mesh, - new_on_plane_nodes, - pointPairPlaneInterception( - *tet4_node_outside_plane, *tet4_nodes_inside_plane[0], plane_point, plane_normal))); - } - auto new_elem_tet4 = std::make_unique(); - new_elem_tet4->set_node(0) = const_cast(new_plane_nodes[0]); - new_elem_tet4->set_node(1) = const_cast(new_plane_nodes[1]); - new_elem_tet4->set_node(2) = const_cast(tet4_nodes_on_plane[0]); - new_elem_tet4->set_node(3) = const_cast(tet4_nodes_inside_plane[0]); - new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id(); - elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4))); - } - else if (tet4_nodes_inside_plane.size() == 2 && tet4_nodes_outside_plane.size() == 1) - { - std::vector new_plane_nodes; - // 2 smaller TET4 elements are created - for (const auto & tet4_node_inside_plane : tet4_nodes_inside_plane) - { - new_plane_nodes.push_back(nonDuplicateNodeCreator( - mesh, - new_on_plane_nodes, - pointPairPlaneInterception( - *tet4_node_inside_plane, *tet4_nodes_outside_plane[0], plane_point, plane_normal))); - } - std::vector new_elems_nodes = {tet4_nodes_inside_plane[0], - tet4_nodes_inside_plane[1], - new_plane_nodes[1], - new_plane_nodes[0], - tet4_nodes_on_plane[0]}; - std::vector> rotated_tet_face_indices; - std::vector> optimized_node_list; - MooseMeshElementConversionUtils::pyramidNodesToTetNodesDeterminer( - new_elems_nodes, rotated_tet_face_indices, optimized_node_list); - - for (unsigned int i = 0; i < optimized_node_list.size(); i++) - { - auto new_elem_tet4 = std::make_unique(); - new_elem_tet4->set_node(0) = const_cast(optimized_node_list[i][0]); - new_elem_tet4->set_node(1) = const_cast(optimized_node_list[i][1]); - new_elem_tet4->set_node(2) = const_cast(optimized_node_list[i][2]); - new_elem_tet4->set_node(3) = const_cast(optimized_node_list[i][3]); - new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id(); - elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4))); - } - } - else - mooseError("Unexpected scenario."); - - mesh.elem_ptr(elem_id)->subdomain_id() = block_id_to_remove; - } - - for (auto & elem_tet4 : elems_tet4) - { - if (new_elements_created) - { - if (elem_tet4->volume() < 0.0) - { - Node * temp = elem_tet4->node_ptr(0); - elem_tet4->set_node(0) = elem_tet4->node_ptr(1); - elem_tet4->set_node(1) = temp; - } - } - // Find the boundary id of the new element - for (unsigned int i = 0; i < 4; i++) - { - const Point side_normal = - (*elem_tet4->side_ptr(i)->node_ptr(1) - *elem_tet4->side_ptr(i)->node_ptr(0)) - .cross(*elem_tet4->side_ptr(i)->node_ptr(2) - *elem_tet4->side_ptr(i)->node_ptr(1)) - .unit(); - for (unsigned int j = 0; j < 4; j++) - { - if (new_elements_created) - { - if (MooseUtils::absoluteFuzzyEqual(side_normal * elem_side_normal_list[j], 1.0)) - { - for (const auto & side_info : elem_side_list[j]) - { - boundary_info.add_side(elem_tet4, i, side_info); - } - } - } - } - if (MooseUtils::absoluteFuzzyEqual(side_normal * _plane_normal, 1.0)) - { - boundary_info.add_side(elem_tet4, i, _cut_face_id); - } - } - } + // Translate the plane point and plane normal to the form of the level set function + _func_level_set = std::make_shared(); + // set FParser internal feature flags + setParserFeatureFlags(_func_level_set); + // The plance is (x - x0) * n_x + (y - y0) * n_y + (z - z0) * n_z = 0 + std::stringstream level_set_ss; + // Let's be conservative about precision here + level_set_ss << std::fixed << std::setprecision(15) << _plane_normal(0) << "*(x-" + << _plane_point(0) << ") + " << _plane_normal(1) << "*(y-" << _plane_point(1) + << ") + " << _plane_normal(2) << "*(z-" << _plane_point(2) << ")"; + + // VERY unlikely to reach this point because we know what we are doing + // But just in case + if (_func_level_set->Parse(level_set_ss.str(), "x,y,z") >= 0) + mooseError("The given plane_point and plane_normal lead to invalid level set.\n", + _func_level_set, + "\nin CutMeshByPlaneGenerator ", + name(), + ".\n", + _func_level_set->ErrorMsg()); + _func_params.resize(3); } diff --git a/test/tests/meshgenerators/cut_mesh_by_level_set_generator/cut_xyzd.i b/test/tests/meshgenerators/cut_mesh_by_level_set_generator/cut_xyzd.i new file mode 100644 index 000000000000..5466674c9620 --- /dev/null +++ b/test/tests/meshgenerators/cut_mesh_by_level_set_generator/cut_xyzd.i @@ -0,0 +1,43 @@ +[Mesh] + [gmg] + type = GeneratedMeshGenerator + dim = 3 + nx = 10 + ny = 10 + nz = 10 + xmin = -1 + ymin = -1 + zmin = -1 + elem_type = HEX8 + [] + [lsc] + type = CutMeshByLevelSetGenerator + input = gmg + level_set = 'x*x+y*y+z*z-0.81' + cut_face_id = 345 + cut_face_name =ls + [] + [xyzd] + type = XYZDelaunayGenerator + boundary = lsc + desired_volume = 1 + [] +[] + +[Executioner] + type = Steady +[] + +[Postprocessors] + [volume] + type = VolumePostprocessor + [] +[] + +[Problem] + solve = false +[] + +[Outputs] + csv = true +[] diff --git a/test/tests/meshgenerators/cut_mesh_by_level_set_generator/gold/cut_xyzd_out.csv b/test/tests/meshgenerators/cut_mesh_by_level_set_generator/gold/cut_xyzd_out.csv new file mode 100644 index 000000000000..007ae58305d1 --- /dev/null +++ b/test/tests/meshgenerators/cut_mesh_by_level_set_generator/gold/cut_xyzd_out.csv @@ -0,0 +1,3 @@ +time,volume +0,0 +1,3.0160040911441 diff --git a/test/tests/meshgenerators/cut_mesh_by_level_set_generator/gold/simple_hex_cut_in.e b/test/tests/meshgenerators/cut_mesh_by_level_set_generator/gold/simple_hex_cut_in.e new file mode 100644 index 0000000000000000000000000000000000000000..6d9142b975b94af10b6d632c909be3e6b6fbc092 GIT binary patch literal 25372 zcmeI3cbt{gwZ8|D4x)k$8}@GOT?f0^3pOkR0?LG{3<_54y(_jPB6f|v*Fh7z#;%EJ zm}rXeCZ^oP#$HN)sqWOL9{U^@n^SIV`t+m(Q?|%1r--83&ZNFpZhPIzx zm<6H7xX}}uCX5(Ax}l+ynVv1p<3~45ZW%pE)%MeU!j$n%qngL7xnwGiu_-V+QpG(2G)bcfy} zCN(!5Gy3Ss%@Zde8@fWLZ8|-35w?2=@gAwS-S0h61@Phz$)qo+2FoYKQ{IYI{y@KcD9zKeh9HMz9Cm#Oq(Z)YUoRcnx5a(xw%C(G~IAOx*ku81>oVN$` zgv{O3=V?>Rv`KAST&d|kbcdn)Y_lO%l{3Cb#c60*`E;BfrJoNt`}vIbGyNs{^`A*S zKbzREPi6kOug9~VAN-zU+QhCaVS7~SCs!VvUc2+=huW^Mw~3+WV|#r5+^y_W2erf8 zuO4mtTJ766O|kUt>}Q;x;cM|9=gomBey4-lVx10ZorCuUeo|L&8|SnQ?#R}+p5OR? zwLX2bj@CId?^(Ii3k`|4x4(VmtRH^b(E7m}OTRJb;)d3j@BQ1p&&+OUoqg?$&3;!d zs(vx()HS~S@kU$THn^ejM|V!x^V2~MnYW!yUFIpRtLypYABMjF=<5xwAD#Qrh&zTi zw7&V-l;>x$o_AVKJpQu1+O44Rx%V%6bqe!7zVYOiS<~CkTlzGnZ~c5($1grP>GDtS zX=r`>(nprufqBdI*H75^%$= z<)aIa|LpVo_TixiAJVw<=lhfED*KfFvDvrkK9#=J^=98npT_L_C;R^K)ce@?Pp)3? zlCkW|C$I4HdQ~@;d}GNsmV9*U$3qV3GmCNQmvw&p*w#(l6`Gc`W(Nm-89jT7S-R-DTa?>noqP zs?PP~`pP-1?sv|0>EG^rXP(k0eH$|_eN&hFQ-9w|U2ENWNWa#ci~94B{vDl*>b}%{ zE;C=o^i3QW!~7tW9JvI-p|r!e!9Nz zzuN!v_g3Bh=X~c}lyzmEcIP?gC7-K&PBJd%Dmr~LPo87Rm;2XH_mk&(-dlC)>U%8j znRfZoH+_=t@O_hWl=tAAwh`6)lXIW6?CeFV**(&*wkShdd7&>V6vM_k2h`&x-){; z`PW+V_2)U;-Z;PKWAgRqGy9(BXIX#FVfB2Lb5#Alk^Rj5$+)~9bM6yQ|Bm*vem^?W z<@ufKj?KR0Tvp@db6eGQ`>ieOsqIy2Fk2%M4O1}O)XWx_0xafBM+umE% z{m;42=ONFL>VD4so~u3&O8>T1b^PD@`sn%*^NZK#8_WFvzqZ|btNOWYTVdv{+I;@L z=yNll$1nQ1FV|K5ewfdF;;~=!dEC+SBYk44>nzWa#AEY3%JpYm>HlS)N9p@zZ8_f^ z{du?D=dFA`^5?6Dx}Ut?^XGxObXBHrKK~uPmy)jzEKcg@_nb|>-T3nUZ!OP{tfyUF zc^-ARj`XQ}PNh%!M>pnsNB#HOIc*i~zW6iCfp2Kh#7y`S%uCN>I4tv0!;4iZFhJ9dP7z+Er{%`;s2>woT7#s}#j`L7B44Pm#jDW*o zB#Z)o1$+b?31grcz7Aty9E=D5Qp!Y_1V_Qq;FG)srob_9EKG%I;NJu|9!`K0;UqX2 zro$<4Dx3zV!x?ZUoCRmYIdCqV2j{~Da3Nd-7sDlRDO?7Z!xbS$& z;ahMW%!KP<7Tf?g!cA~9+yb}4ZE!o>0pEr@;V!rv?ty#ZKDZwqfCu3rco-gmN8vGe z9G-wD;VEc^+3++x1K)w~!n5!kJP$9xi|`V>46nee@I819UWYf}O?V65hIimycn>;3 zXXpY86y3GNqPwAuMR&(H7Tp8iShT&2sqNhpdO`1^`;=I;y^TdX&RBHT7HuzMYUk1K z3kw$Qcw^C7TXa8i#-bO(Hx_L#V`}?2-@>p+(aviu+WCw{JI+{i))s9qW9p^Bd$0^F zTcejNvE&D#jYV6#JUG_aAW*yCE5M3i&c6D_qV*lGHs_d?Amh|4mss+vlvwntXk%LI zji;Y;8dq!8u3ha}SPfQ(%%ku8YOUHDwSG@n12Rs%W{D-gR*6NgjW(vW-gx?XZN}AF zwb!DyHoey=zA<~LzXqPUbz!~YYmIw-?yJ7BX#MqzPHg!S+W_mB6-r*~I*ePZRd*@c z-qzRV&AF!3>t{XYHU#Sv*Ebfe@Aat7Ic6isIQ7OQrd3;)b-HHbYOOl=bCcq)hBdEs zKaH!kYClgI@7UVBIp<7$)}uDJDOjJlzUx9kUr^oO<&T)2gk@I$g7IwN|}G z(XpP5Eg-$NEO}%4sne>{FZHg+_*USvXKQFIzSejT@IL7qi`E}pbYh-4W7~jZa=*3i zgK=xM>b^zW+xptPIoFhW{jA5_wqSkY`o^O5y&knW$7}}~r{2EAv})_JPSYYnWtF|ud zbj`-qTD9jv9qZW`0_n9&$s5y8omQQGsdqibb02oa-wnJcIcsXKVRy*=GpBDXwfcJ$ zotXWMS!)mX#X9HT6S6z0T8!JK0b0OQ8=2Z36v&itvjpYa31HDrFRHO3vURd+3V zRjl>3d2^jh-nhPXjyHD@)W-E!0qeACYt;IV83q}rKDflRYU{F2*KAy^Rl9%cL%=!G zU+Z||&ZkxPDSBvYY= z)mpXtlXVZrPk-mrH|Bh5{Slx}%=d#aYn{*a3IsMdX z)z+9->yIkBd-3gKUTdzlUOgJjdk^(9UafB*^X81@{A*o%#%pV@#p}#^5V_vO4uIPTnw~u+Px!QX5crf3+_!+O(w~u-2jCoJh zTF*kpYt>m_FMRi42uuLAdtgl68z#aSNKQYsTD3Li)%uf)UZnW;F|RdOTdzI}%r9L0 zj92U1$GmmMb}3rx+Kne)do2qQ&;B2cuimZ1qIWH^=*ehf>VD7yQ^2w2^i!)WuG;e+=vn*6&e#W6}D@7OhoJEqcM?+sAsXxs117Jq^seHvNoO>)Xe? zIb*qAt=F9KT6Nar8273+<~10z#(Q-f>zt@?zb7b(7dtk;^$cp&mRcAfc>|Pi18jM-vdOU-B!%0OSRbp!AJsBo}wdV9w zt5sWLUR!(3xz^ge*Pii?ac^s5dzY9s-JmNh0R4*|P-4-GqK!o_hHos|=ZUdsdl^&P zyASk+#fx5|#G>tOEZT9#qO-PWdl^$ZkG}V8$)X)^EIMn89!SnuwCCAaw7ra}?c;pz z`BFtYud!(7GZyVQW6@b#w7ra}-N)&03V1K}EB?@;wdST5t?!sqVLC*gR$`8~&jD~c zSZmE0VBDB~>a=QW%%|Rd#?J)TFdfWmtue0Fs=Z(8gTeaRyt%HBdVTBE=FS4^6W1SB zv{r46d40#64H>6Cr^K{s>$1-2_{P=R>EQmU&jsg5f34$*Gd zc(uNL%$qZIM$y{SA>*~`tk3Ho4oAY}pmuMDgWA2g0%+sC}tTy4GjS}?y@@iShnKM>4YXUu!5)_N8)UaQXfvKJGHeFL(8?u*)e z`zE-D$?2z7tG33xTK`)`FI9Z|nAe)Ctyf%FYat7n4w1>r`RTJ*Hy z-&C~LoHf>~^=~eE$>Q55{mj+Yt8XcJ*QTHGYJK}S-kh;quhwhMc&$3?bBudc8}k~B zS>wIB6;1%xcw+HyD_U#L8uM!X+l%(T>)Xe?R$W`Kz5~pA9r_us*0+y&bH==8wbs4L zc&$3?vnKnWnAc!TowIl{{Dy(pXpr@~#uH+BZx4au8RTW3st4-9~F;9f9pOh0v6wKe8be|Cv$jahRaECCn5 z{b1afe(JPpYs{zq{1Vq1v*v-~8@sS**N|&_5MO;+iK%@yJ_MISa{8&&s;x1vb zYafBstMyaw827HW#(gzr-J|drJPy`+jrzu-_4^lX&YJAM)_aigTD5&0FS&P#xfWw; zujL8I^_bH)=2*3U#%k4BpX<`MkM&x0#y^Sw6tqI_pT4nZ{l$wmw*+`U&1*dq8Lw5_ z*SzE`CFUL&Q@g&|;F_#8r#~B_^)pth&iY)JzJ09Ms%L}RwO$BUz|%!PP-4+ugRy9@ z!I;|q&wxGcWlrB%w0_3wJ7z}7Yn?mw+S;{vKXQJ*gRj1}#MA@fyKoI8r=MD_+8Xm( z$2?oK_BlwsT0ixUc^GPI9x5^Go`)CUMX=WW)i)Nczhu$ota+?xt!K?X)~of?&zSSq z*0>g9)_E;2L9WM~zA?wD_0wOg&iY)JzJ09Msx$s&{8!*raNXWJePhx3OBHR-nw+_u zHTzhv)=$6e#q1JSd+qLnb*}Gw;F^-tHx{j*{#td`=eqRmW4%_L@ve0Sd=p*+wP(;+ zbk1|MYd5C0|Lb7S+PyIrt)H>_j`>!}Yn?mw+S;{vKdu4K;~Su!1;(OpD6!}_(Z-@Z z55}VHWlZgUyahLdy>BW0+eK^5Sz})7n4FoMmDH>CQ}39Z$JDqMW6`eFShQ;}rgpF1 zf!o1(?KLtdO%O;1-+pU^o0eXA1nk5!y?ch2Ed}Q7%UDx zpVx!+VFTC@HiC^|6WA0sgUw+J*b=sat)USH!#1!jYzN!J4zMHa1UthJ*adcl-C%dv z1NMZyU~kw5_JyIaAM6hYz=3cO41%905nd7-)vC!&n#x z<6#0!gh_A|91WA91*X6;a4bxPX>c4I4=2Eha1xvh)8Q026;6ZG;S4wv&VsYy95@%w zgY)46xDYOai{TQu6fT3y;R=`mSHe|rHCzMN!Z*NYgJr0k9}628+WIup|tGrC@1T z29|~8U=S=1E5M4d608iXz^bqstPX3yny?nE4eP+yU|m=b)`tyXL)Zv5hD~5ou9NvI z&wb|q-}t-Hoi<~h&0!1J61IY^p%DhdHn1&h2iwCAup{gQJHrsz1$KqqV0YL9_JqA) zZ`cR+g`uz?><Q zcnMyHSKw9n9=ryx!yE7>yajK=JMb>N2j7Pu!29q+_!0aVegZ#*Iq(7e41NwD!bk86 z_$B-behnYPZ{QR75BM$o4nBq7!yn*}@F)0B_%r+kK7+r)-{8OC@9+=!C;WHYYrk_R zY-i{K3qV)s2Hl|t^n_l}8~Q+BSP=TbLa;C_0{vkCEDDRk;;;lP2?Jp%SQ?grWnnoO z1k1w;up+DkE5j01_Reu{r<$YRlhrN zP3^uXaShe)NUVRZ??;EA4~4^^2@Z!N;7Aw)%`g_m!FZSiEieU+fn#AB90$k4iEt8} z4AbEhI1NsRGvG`(8_t1q;XF7WE`W>RVz>k@h0EXym;qP9Rd6j_2Q%S%m<2b&O>i^Z z26w_;a5vlo_rd+}AUp&Q!z1u0JO)p|lkgNg4bQ@J@I1T#FTu<33cLny!Q1c-ybIrl zAHe(YWB34m20w=n;TP~r_!ayHeg~hz@8J*dNB9%`8U6x)g}=ey;UDl%KG&U~D|CbI z&;xowZ|DR4pg#pjHzslJ0A!{@xr z^L@klsl01mz`h8apLxD#-XQ)aaGomfn0K-70q3dR`{gIZehMGJ&w#U4eXsl)|6|~6 zRo*4P#d5|f?-9;cyLSfXD({PS?+MOVLbW7q^XgUx~SInU?&w#2uC9biY;33i4dunX)4 zy8~zSOFhrK@Oi7HGkSqSZCpR`*pBk{NG`RteG{v%Nj$bbuTCy=)av-L_2Uzdt