LCOV - code coverage report
Current view: top level - src/meshgenerators - Boundary2DDelaunayGenerator.C (source / functions) Hit Total Coverage
Test: idaholab/moose framework: 909fe5 Lines: 228 241 94.6 %
Date: 2025-08-29 20:01:24 Functions: 9 9 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : //* This file is part of the MOOSE framework
       2             : //* https://www.mooseframework.org
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "Boundary2DDelaunayGenerator.h"
      11             : 
      12             : #include "MooseMeshUtils.h"
      13             : #include "CastUniquePointer.h"
      14             : 
      15             : #include "libmesh/boundary_info.h"
      16             : #include "libmesh/poly2tri_triangulator.h"
      17             : #include "libmesh/mesh_triangle_holes.h"
      18             : #include "libmesh/mesh_modification.h"
      19             : #include "libmesh/parallel_algebra.h"
      20             : 
      21             : registerMooseObject("MooseApp", Boundary2DDelaunayGenerator);
      22             : 
      23             : InputParameters
      24       14675 : Boundary2DDelaunayGenerator::validParams()
      25             : {
      26       14675 :   InputParameters params = SurfaceDelaunayGeneratorBase::validParams();
      27       14675 :   params += FunctionParserUtils<false>::validParams();
      28             : 
      29       29350 :   params.addClassDescription(
      30             :       "Mesh generator that convert a 2D surface given as one or a few boundaries of a 3D mesh into "
      31             :       "a 2D mesh using Delaunay triangulation.");
      32       58700 :   params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
      33       58700 :   params.addRequiredParam<std::vector<BoundaryName>>("boundary_names", "The boundaries to be used");
      34       44025 :   params.addParam<std::vector<std::vector<BoundaryName>>>(
      35             :       "hole_boundary_names",
      36       29350 :       std::vector<std::vector<BoundaryName>>(),
      37             :       "The optional boundaries to be used as the holes in the mesh during triangulation. Note that "
      38             :       "this is a vector of vectors, which allows each hole to be defined as a combination of "
      39             :       "multiple boundaries.");
      40             : 
      41       58700 :   params.addParam<std::string>(
      42             :       "level_set",
      43             :       "Level set used to achieve more accurate reverse projection compared to interpolation.");
      44       44025 :   params.addParam<unsigned int>(
      45             :       "max_level_set_correction_iterations",
      46       29350 :       3,
      47             :       "Maximum number of iterations to correct the nodes based on the level set function.");
      48       73375 :   params.addRangeCheckedParam<Real>(
      49             :       "max_angle_deviation",
      50       29350 :       60.0,
      51             :       "max_angle_deviation>0 & max_angle_deviation<90",
      52             :       "Maximum angle deviation from the global average normal vector in the input mesh.");
      53             : 
      54       44025 :   params.addParam<BoundaryName>(
      55             :       "output_external_boundary_name",
      56             :       "",
      57             :       "The optional name of the external boundary of the mesh to generate. If not provided, the "
      58             :       "external boundary will be have a trivial id of 0.");
      59             : 
      60       14675 :   return params;
      61           0 : }
      62             : 
      63         205 : Boundary2DDelaunayGenerator::Boundary2DDelaunayGenerator(const InputParameters & parameters)
      64             :   : SurfaceDelaunayGeneratorBase(parameters),
      65             :     FunctionParserUtils<false>(parameters),
      66         205 :     _input(getMesh("input")),
      67         410 :     _boundary_names(getParam<std::vector<BoundaryName>>("boundary_names")),
      68         410 :     _hole_boundary_names(getParam<std::vector<std::vector<BoundaryName>>>("hole_boundary_names")),
      69         205 :     _max_level_set_correction_iterations(
      70         410 :         getParam<unsigned int>("max_level_set_correction_iterations")),
      71         410 :     _max_angle_deviation(getParam<Real>("max_angle_deviation")),
      72         820 :     _output_external_boundary_name(getParam<BoundaryName>("output_external_boundary_name"))
      73             : {
      74         615 :   if (isParamValid("level_set"))
      75             :   {
      76          49 :     _func_level_set = std::make_shared<SymFunction>();
      77             :     // set FParser internal feature flags
      78          49 :     setParserFeatureFlags(_func_level_set);
      79         147 :     if (isParamValid("constant_names") && isParamValid("constant_expressions"))
      80           0 :       addFParserConstants(_func_level_set,
      81             :                           getParam<std::vector<std::string>>("constant_names"),
      82             :                           getParam<std::vector<std::string>>("constant_expressions"));
      83         245 :     if (_func_level_set->Parse(getParam<std::string>("level_set"), "x,y,z") >= 0)
      84           0 :       mooseError("Invalid function f(x,y,z)\n",
      85           0 :                  _func_level_set,
      86             :                  "\nin CutMeshByLevelSetGenerator ",
      87           0 :                  name(),
      88             :                  ".\n",
      89           0 :                  _func_level_set->ErrorMsg());
      90             : 
      91          49 :     _func_params.resize(3);
      92             :   }
      93         205 : }
      94             : 
      95             : std::unique_ptr<MeshBase>
      96         115 : Boundary2DDelaunayGenerator::generate()
      97             : {
      98         115 :   std::unique_ptr<MeshBase> mesh_3d = std::move(_input);
      99             : 
     100             :   // Generate a new 2D block based on the sidesets
     101         115 :   const auto new_block_id = MooseMeshUtils::getNextFreeSubdomainID(*mesh_3d);
     102             :   try
     103             :   {
     104         242 :     MooseMeshUtils::createSubdomainFromSidesets(
     105         349 :         mesh_3d, _boundary_names, new_block_id, SubdomainName(), type());
     106             :   }
     107           4 :   catch (MooseException & e)
     108             :   {
     109           8 :     if (((std::string)e.what()).compare(0, 12, "The sideset ") == 0)
     110           8 :       paramError("boundary_names", e.what());
     111             :     else
     112           0 :       mooseError(e.what());
     113           0 :   }
     114             : 
     115             :   // If holes are provided, we need to create new blocks for them too
     116         111 :   std::vector<subdomain_id_type> hole_block_ids;
     117         128 :   for (const auto & hole : _hole_boundary_names)
     118             :   {
     119          21 :     hole_block_ids.push_back(MooseMeshUtils::getNextFreeSubdomainID(*mesh_3d));
     120             :     try
     121             :     {
     122          71 :       MooseMeshUtils::createSubdomainFromSidesets(
     123          67 :           mesh_3d, hole, hole_block_ids.back(), SubdomainName(), type());
     124             :     }
     125           4 :     catch (MooseException & e)
     126             :     {
     127           8 :       if (((std::string)e.what()).compare(0, 12, "The sideset ") == 0)
     128           8 :         paramError("hole_boundary_names", e.what());
     129             :       else
     130           0 :         mooseError(e.what());
     131           0 :     }
     132             :   }
     133             : 
     134             :   // Create a 2D mesh form the 2D block
     135         107 :   auto mesh_2d = buildMeshBaseObject();
     136         321 :   MooseMeshUtils::convertBlockToMesh(mesh_3d, mesh_2d, {std::to_string(new_block_id)});
     137             :   // If holes are provided, we need to create a 2D mesh for each hole
     138         107 :   std::vector<std::unique_ptr<MeshBase>> hole_meshes_2d;
     139         124 :   for (const auto & hole_block_id : hole_block_ids)
     140             :   {
     141          17 :     hole_meshes_2d.push_back(buildMeshBaseObject());
     142          68 :     MooseMeshUtils::convertBlockToMesh(
     143          34 :         mesh_3d, hole_meshes_2d.back(), {std::to_string(hole_block_id)});
     144             :   }
     145             : 
     146             :   // We do not need the 3D mesh anymore
     147         107 :   mesh_3d->clear();
     148             : 
     149         206 :   return General2DDelaunay(mesh_2d, hole_meshes_2d);
     150         223 : }
     151             : 
     152             : Point
     153        8029 : Boundary2DDelaunayGenerator::elemNormal(const Elem & elem)
     154             : {
     155             :   mooseAssert(elem.n_vertices() == 3 || elem.n_vertices() == 4, "unsupported element type.");
     156             :   // Only the first three vertices are used to calculate the normal vector
     157        8029 :   const Point & p0 = *elem.node_ptr(0);
     158        8029 :   const Point & p1 = *elem.node_ptr(1);
     159        8029 :   const Point & p2 = *elem.node_ptr(2);
     160             : 
     161        8029 :   if (elem.n_vertices() == 4)
     162             :   {
     163           0 :     const Point & p3 = *elem.node_ptr(3);
     164           0 :     return ((p2 - p0).cross(p3 - p1)).unit();
     165             :   }
     166             : 
     167        8029 :   return ((p2 - p1).cross(p0 - p1)).unit();
     168             : }
     169             : 
     170             : Point
     171         103 : Boundary2DDelaunayGenerator::meshNormal2D(const MeshBase & mesh)
     172             : {
     173         103 :   Point mesh_norm = Point(0.0, 0.0, 0.0);
     174         103 :   Real mesh_area = 0.0;
     175             : 
     176             :   // Check all the elements' normal vectors
     177        2787 :   for (const auto & elem : mesh.active_local_element_ptr_range())
     178             :   {
     179        2684 :     const Real elem_area = elem->volume();
     180        2684 :     mesh_norm += elemNormal(*elem) * elem_area;
     181        2684 :     mesh_area += elem_area;
     182         103 :   }
     183         103 :   mesh.comm().sum(mesh_norm);
     184         103 :   mesh.comm().sum(mesh_area);
     185         103 :   mesh_norm /= mesh_area;
     186         103 :   return mesh_norm.unit();
     187             : }
     188             : 
     189             : Real
     190         103 : Boundary2DDelaunayGenerator::meshNormalDeviation2D(const MeshBase & mesh, const Point & global_norm)
     191             : {
     192         103 :   Real max_deviation(0.0);
     193             :   // Check all the elements' deviation from the global normal vector
     194        2787 :   for (const auto & elem : mesh.active_local_element_ptr_range())
     195             :   {
     196        2684 :     const Real elem_deviation = std::acos(global_norm * elemNormal(*elem)) / M_PI * 180.0;
     197        2684 :     max_deviation = std::max(max_deviation, elem_deviation);
     198         103 :   }
     199         103 :   mesh.comm().max(max_deviation);
     200         103 :   return max_deviation;
     201             : }
     202             : 
     203             : Real
     204       12139 : Boundary2DDelaunayGenerator::levelSetEvaluator(const Point & point)
     205             : {
     206       48556 :   return evaluate(_func_level_set, std::vector<Real>({point(0), point(1), point(2), 0}));
     207             : }
     208             : 
     209             : void
     210        1320 : Boundary2DDelaunayGenerator::levelSetCorrection(Node & node)
     211             : {
     212             :   // Based on the given level set, we try to move the node in its normal direction
     213        1320 :   const Real diff = libMesh::TOLERANCE * 10.0; // A small value to perturb the node
     214        1320 :   const Real original_eval = levelSetEvaluator(node);
     215        1320 :   const Real xp_eval = levelSetEvaluator(node + Point(diff, 0.0, 0.0));
     216        1320 :   const Real yp_eval = levelSetEvaluator(node + Point(0.0, diff, 0.0));
     217        1320 :   const Real zp_eval = levelSetEvaluator(node + Point(0.0, 0.0, diff));
     218        1320 :   const Real xm_eval = levelSetEvaluator(node - Point(diff, 0.0, 0.0));
     219        1320 :   const Real ym_eval = levelSetEvaluator(node - Point(0.0, diff, 0.0));
     220        1320 :   const Real zm_eval = levelSetEvaluator(node - Point(0.0, 0.0, diff));
     221        1320 :   const Point grad = Point((xp_eval - xm_eval) / (2.0 * diff),
     222        1320 :                            (yp_eval - ym_eval) / (2.0 * diff),
     223        1320 :                            (zp_eval - zm_eval) / (2.0 * diff));
     224        1320 :   const Real xyz_diff = -original_eval / grad.contract(grad);
     225        1320 :   node(0) += xyz_diff * grad(0);
     226        1320 :   node(1) += xyz_diff * grad(1);
     227        1320 :   node(2) += xyz_diff * grad(2);
     228        1320 : }
     229             : 
     230             : std::unique_ptr<MeshBase>
     231         107 : Boundary2DDelaunayGenerator::General2DDelaunay(
     232             :     std::unique_ptr<MeshBase> & mesh_2d, std::vector<std::unique_ptr<MeshBase>> & hole_meshes_2d)
     233             : {
     234             :   // If a level set is provided, we need to check if the nodes in the original 2D mesh match the
     235             :   // level set
     236         107 :   if (_func_level_set)
     237             :   {
     238         949 :     for (const auto & node : mesh_2d->node_ptr_range())
     239             :     {
     240         904 :       if (std::abs(levelSetEvaluator(*node)) > libMesh::TOLERANCE)
     241             :       {
     242           8 :         paramError("level_set",
     243             :                    "The level set function does not match the nodes in the given boundary of the "
     244             :                    "input mesh.");
     245             :       }
     246          45 :     }
     247             :   }
     248             : 
     249             :   // Easier to work with a TRI3 mesh
     250             :   // all_tri() also prepares the mesh for use
     251         103 :   mesh_2d->prepare_for_use();
     252         103 :   MeshTools::Modification::all_tri(*mesh_2d);
     253         103 :   const auto mesh_2d_ext_bdry = MooseMeshUtils::getNextFreeBoundaryID(*mesh_2d);
     254        3499 :   for (const auto & elem : mesh_2d->active_element_ptr_range())
     255       13584 :     for (const auto & i_side : elem->side_index_range())
     256       10188 :       if (elem->neighbor_ptr(i_side) == nullptr)
     257        1675 :         mesh_2d->get_boundary_info().add_side(elem, i_side, mesh_2d_ext_bdry);
     258             : 
     259             :   // Create a clone of the 2D mesh to be used for the 1D mesh generation
     260         103 :   auto mesh_2d_dummy = dynamic_pointer_cast<MeshBase>(mesh_2d->clone());
     261             :   // Generate a new 1D block based on the external boundary
     262         103 :   const auto new_block_id_1d = MooseMeshUtils::getNextFreeSubdomainID(*mesh_2d_dummy);
     263             : 
     264         412 :   MooseMeshUtils::createSubdomainFromSidesets(
     265         309 :       mesh_2d_dummy, {std::to_string(mesh_2d_ext_bdry)}, new_block_id_1d, SubdomainName(), type());
     266             : 
     267             :   // Create a 1D mesh form the 1D block
     268         103 :   auto mesh_1d = buildMeshBaseObject();
     269         309 :   MooseMeshUtils::convertBlockToMesh(mesh_2d_dummy, mesh_1d, {std::to_string(new_block_id_1d)});
     270         103 :   mesh_2d_dummy->clear();
     271             : 
     272             :   // If we have holes, we need to create a 1D mesh for each hole
     273         103 :   std::vector<std::unique_ptr<MeshBase>> hole_meshes_1d;
     274         116 :   for (auto & hole_mesh_2d : hole_meshes_2d)
     275             :   {
     276             :     // As we do not need these holes for reverse projection, we do not need to convert them to TRI3
     277             :     // meshes, but we still need to create a 1D mesh for each hole
     278          13 :     hole_mesh_2d->find_neighbors();
     279          13 :     const auto hole_mesh_2d_ext_bdry = MooseMeshUtils::getNextFreeBoundaryID(*hole_mesh_2d);
     280          65 :     for (const auto & elem : hole_mesh_2d->active_element_ptr_range())
     281         260 :       for (const auto & i_side : elem->side_index_range())
     282         208 :         if (elem->neighbor_ptr(i_side) == nullptr)
     283         117 :           hole_mesh_2d->get_boundary_info().add_side(elem, i_side, mesh_2d_ext_bdry);
     284          13 :     const auto new_hole_block_id_1d = MooseMeshUtils::getNextFreeSubdomainID(*hole_mesh_2d);
     285          52 :     MooseMeshUtils::createSubdomainFromSidesets(hole_mesh_2d,
     286          13 :                                                 {std::to_string(hole_mesh_2d_ext_bdry)},
     287             :                                                 new_hole_block_id_1d,
     288          26 :                                                 SubdomainName(),
     289          13 :                                                 type());
     290             :     // Create a 1D mesh form the 1D block
     291          13 :     hole_meshes_1d.push_back(buildMeshBaseObject());
     292          52 :     MooseMeshUtils::convertBlockToMesh(
     293          26 :         hole_mesh_2d, hole_meshes_1d.back(), {std::to_string(new_hole_block_id_1d)});
     294          13 :     hole_mesh_2d->clear();
     295             :   }
     296             : 
     297             :   // Find centroid of the 2D mesh
     298         103 :   const Point centroid = MooseMeshUtils::meshCentroidCalculator(*mesh_2d);
     299             :   // calculate an average normal vector of the 2D mesh
     300         103 :   const Point mesh_norm = meshNormal2D(*mesh_2d);
     301             :   // Check the deviation of the mesh normal vector from the global average normal vector
     302         103 :   if (meshNormalDeviation2D(*mesh_2d, mesh_norm) > _max_angle_deviation)
     303           4 :     paramError("boundary_names",
     304             :                "The normal vector of some elements in the 2D mesh deviates too much from the "
     305           4 :                "global average normal vector. The maximum deviation is " +
     306           8 :                    std::to_string(_max_angle_deviation) +
     307             :                    ". Consider dividing the boundary into several parts to "
     308             :                    "reduce the angle deviation.");
     309             : 
     310             :   // Move both 2d and 1d meshes to the centroid of the 2D mesh
     311          99 :   MeshTools::Modification::translate(*mesh_1d, -centroid(0), -centroid(1), -centroid(2));
     312          99 :   MeshTools::Modification::translate(*mesh_2d, -centroid(0), -centroid(1), -centroid(2));
     313             :   // Alsop need to translate the 1D hole meshes if applicable
     314         108 :   for (auto & hole_mesh_1d : hole_meshes_1d)
     315           9 :     MeshTools::Modification::translate(*hole_mesh_1d, -centroid(0), -centroid(1), -centroid(2));
     316             : 
     317             :   // Calculate the Euler angles to rotate the meshes so that the 2D mesh is close to the XY plane
     318             :   // (i.e., the normal vector of the 2D mesh is aligned with the Z axis)
     319          99 :   const Real theta = std::acos(mesh_norm(2)) / M_PI * 180.0;
     320             :   const Real phi =
     321          99 :       (MooseUtils::absoluteFuzzyLessThan(mesh_norm(2), 1.0) ? std::atan2(mesh_norm(1), mesh_norm(0))
     322          99 :                                                             : 0.0) /
     323          99 :       M_PI * 180.0;
     324          99 :   MeshTools::Modification::rotate(*mesh_1d, 90.0 - phi, theta, 0.0);
     325          99 :   MeshTools::Modification::rotate(*mesh_2d, 90.0 - phi, theta, 0.0);
     326             :   // Also rotate the 1D hole meshes if applicable
     327         108 :   for (auto & hole_mesh_1d : hole_meshes_1d)
     328           9 :     MeshTools::Modification::rotate(*hole_mesh_1d, 90.0 - phi, theta, 0.0);
     329             : 
     330             :   // Clone the 2D mesh to be used for reverse projection later
     331          99 :   auto mesh_2d_xyz = dynamic_pointer_cast<MeshBase>(mesh_2d->clone());
     332             : 
     333             :   // Project the 2D mesh to the XY plane so that XYDelaunay can be used
     334        2538 :   for (const auto & node : mesh_2d->node_ptr_range())
     335        2538 :     (*node)(2) = 0;
     336             :   // Project the 1D mesh to the XY plane as well
     337        1575 :   for (const auto & node : mesh_1d->node_ptr_range())
     338        1575 :     (*node)(2) = 0;
     339             :   // Also project the 1D hole meshes to the XY plane if applicable
     340         108 :   for (auto & hole_mesh_1d : hole_meshes_1d)
     341             :   {
     342          81 :     for (const auto & node : hole_mesh_1d->node_ptr_range())
     343          81 :       (*node)(2) = 0;
     344             :   }
     345             : 
     346          99 :   std::vector<libMesh::TriangulatorInterface::MeshedHole> meshed_holes;
     347          99 :   std::vector<libMesh::TriangulatorInterface::Hole *> triangulator_hole_ptrs(hole_meshes_1d.size());
     348          99 :   meshed_holes.reserve(hole_meshes_1d.size());
     349         108 :   for (auto hole_i : index_range(hole_meshes_1d))
     350             :   {
     351           9 :     hole_meshes_1d[hole_i]->prepare_for_use();
     352           9 :     meshed_holes.emplace_back(*hole_meshes_1d[hole_i]);
     353           9 :     triangulator_hole_ptrs[hole_i] = &meshed_holes.back();
     354             :   }
     355             : 
     356             :   // Finally, triangulation
     357             :   std::unique_ptr<UnstructuredMesh> mesh =
     358          99 :       dynamic_pointer_cast<UnstructuredMesh>(std::move(mesh_1d));
     359             : 
     360          99 :   Poly2TriTriangulator poly2tri(*mesh);
     361          99 :   poly2tri.triangulation_type() = TriangulatorInterface::PSLG;
     362             : 
     363          99 :   poly2tri.set_interpolate_boundary_points(0);
     364          99 :   poly2tri.set_refine_boundary_allowed(false);
     365          99 :   poly2tri.set_verify_hole_boundaries(false);
     366          99 :   poly2tri.desired_area() = 0;
     367          99 :   poly2tri.minimum_angle() = 0; // Not yet supported
     368          99 :   poly2tri.smooth_after_generating() = false;
     369          99 :   if (!triangulator_hole_ptrs.empty())
     370           9 :     poly2tri.attach_hole_list(&triangulator_hole_ptrs);
     371             :   // Future TODO: correct the area function based on the local normal vector
     372          99 :   if (_use_auto_area_func)
     373          90 :     poly2tri.set_auto_area_function(this->comm(),
     374          90 :                                     _auto_area_function_num_points,
     375          90 :                                     _auto_area_function_power,
     376          90 :                                     _auto_area_func_default_size,
     377          90 :                                     _auto_area_func_default_size_dist);
     378          99 :   poly2tri.triangulate();
     379             : 
     380             :   // Reverse the projection based on the original 2D mesh
     381        2760 :   for (const auto & node : mesh->node_ptr_range())
     382             :   {
     383        2661 :     bool node_mod = false;
     384             :     // Try to find the element in mesh_2d that contains the new node
     385       47083 :     for (const auto & elem : mesh_2d->active_element_ptr_range())
     386             :     {
     387       47083 :       if (elem->contains_point(Point((*node)(0), (*node)(1), 0.0)))
     388             :       {
     389             :         // Element id
     390        2661 :         const auto elem_id = elem->id();
     391             :         // element in xyz_in_xyz
     392        2661 :         const Elem & elem_xyz = *mesh_2d_xyz->elem_ptr(elem_id);
     393             : 
     394        2661 :         const Point elem_normal = elemNormal(elem_xyz);
     395        2661 :         const Point & elem_p = *mesh_2d_xyz->elem_ptr(elem_id)->node_ptr(0);
     396             : 
     397             :         // if the x and y values of the node is the same as the elem_p's first node, we can just
     398             :         // move it to that node's position
     399        3561 :         if (MooseUtils::absoluteFuzzyEqual((*node)(0), elem_p(0)) &&
     400         900 :             MooseUtils::absoluteFuzzyEqual((*node)(1), elem_p(1)))
     401             :         {
     402         189 :           (*node)(2) = elem_p(2);
     403         189 :           node_mod = true;
     404         189 :           break;
     405             :         }
     406             :         // Otherwise, we need to find a position inside the 2D element
     407             :         // It has the same x and y coordinates as the node in the projected mesh;
     408        2472 :         (*node)(2) = elem_p(2) - (((*node)(0) - elem_p(0)) * elem_normal(0) +
     409        2472 :                                   ((*node)(1) - elem_p(1)) * elem_normal(1)) /
     410        2472 :                                      elem_normal(2);
     411        2472 :         node_mod = true;
     412        2472 :         break;
     413             :       }
     414        2661 :     }
     415        2661 :     if (!node_mod)
     416           0 :       mooseError("Node not found in mesh_in_xy");
     417          99 :   }
     418             : 
     419             :   // Rotate the mesh back
     420          99 :   MeshTools::Modification::rotate(*mesh, 0.0, -theta, phi - 90.0);
     421             :   // Translate the mesh back
     422          99 :   MeshTools::Modification::translate(*mesh, centroid(0), centroid(1), centroid(2));
     423             : 
     424             :   // Correct the nodes based on the level set function
     425          99 :   if (_func_level_set)
     426             :   {
     427        1160 :     for (const auto & node : mesh->node_ptr_range())
     428             :     {
     429        1115 :       unsigned int iter_ct = 0;
     430        4430 :       while (iter_ct < _max_level_set_correction_iterations &&
     431        1995 :              std::abs(levelSetEvaluator(*node)) > libMesh::TOLERANCE * libMesh::TOLERANCE)
     432             :       {
     433        1320 :         levelSetCorrection(*node);
     434        1320 :         ++iter_ct;
     435             :       }
     436          45 :     }
     437             :   }
     438             : 
     439             :   // Assign the external boundary name if provided
     440          99 :   if (!_output_external_boundary_name.empty())
     441             :   {
     442             :     const std::vector<BoundaryID> output_boundary_id =
     443          54 :         MooseMeshUtils::getBoundaryIDs(*mesh, {_output_external_boundary_name}, true);
     444             : 
     445          18 :     libMesh::MeshTools::Modification::change_boundary_id(*mesh, 0, output_boundary_id[0]);
     446          18 :     mesh->get_boundary_info().sideset_name(output_boundary_id[0]) = _output_external_boundary_name;
     447          18 :   }
     448             : 
     449          99 :   mesh->set_isnt_prepared();
     450             : 
     451         198 :   return mesh;
     452         349 : }

Generated by: LCOV version 1.14