Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       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 "ProjectSideSetOntoLevelSetGenerator.h"
      11             : #include "MooseMeshUtils.h"
      12             : #include "CastUniquePointer.h"
      13             : 
      14             : #include "libmesh/elem.h"
      15             : #include "libmesh/boundary_info.h"
      16             : #include "libmesh/mesh_base.h"
      17             : 
      18             : // C++ includes
      19             : #include <cmath>
      20             : 
      21             : registerMooseObject("MooseApp", ProjectSideSetOntoLevelSetGenerator);
      22             : 
      23             : InputParameters
      24        3077 : ProjectSideSetOntoLevelSetGenerator::validParams()
      25             : {
      26        3077 :   InputParameters params = MeshGenerator::validParams();
      27        3077 :   params += FunctionParserUtils<false>::validParams();
      28             : 
      29       12308 :   params.addRequiredParam<Point>("direction", "Projection direction");
      30       12308 :   params.addParam<SubdomainName>(
      31             :       "subdomain_name", "projected", "Name of the subdomain of the projected surface");
      32        9231 :   params.addParam<Real>("max_search_distance",
      33        6154 :                         1e6,
      34             :                         "Maximum distance from the sideset to search for the projected point in "
      35             :                         "the projection direction.");
      36             : 
      37             :   // Source domain parameters
      38       12308 :   params.addRequiredParam<MeshGeneratorName>(
      39             :       "input", "The input mesh generator in which the sideset exists.");
      40       12308 :   params.addRequiredParam<BoundaryName>("sideset", "The name of the sideset to project");
      41             : 
      42             :   // Surface definition
      43       12308 :   params.addRequiredParam<std::string>(
      44             :       "level_set",
      45             :       "Level set to define the surface to project onto as a function of x, y, and z. The surface "
      46             :       "should be 'in front of' the sideset when following the projection direction.");
      47       12308 :   params.addParam<std::vector<std::string>>(
      48             :       "constant_names", {}, "Vector of constants used in the parsed function");
      49       12308 :   params.addParam<std::vector<std::string>>(
      50             :       "constant_expressions",
      51             :       {},
      52             :       "Vector of values for the constants in constant_names (can be an FParser expression)");
      53       12308 :   params.addParamNamesToGroup("level_set constant_names constant_expressions", "Level set surface");
      54             : 
      55        3077 :   params.addClassDescription(
      56             :       "Projects a sideset onto a surface defined by a level set and creates a surface mesh.");
      57             : 
      58        3077 :   return params;
      59           0 : }
      60             : 
      61           8 : ProjectSideSetOntoLevelSetGenerator::ProjectSideSetOntoLevelSetGenerator(
      62           8 :     const InputParameters & parameters)
      63             :   : MeshGenerator(parameters),
      64             :     FunctionParserUtils<false>(parameters),
      65           8 :     _input_name(getParam<MeshGeneratorName>("input")),
      66           8 :     _input(getMeshByName(_input_name)),
      67          16 :     _level_set(getParam<std::string>("level_set")),
      68          40 :     _max_search_distance(getParam<Real>("max_search_distance"))
      69             : {
      70             :   // Create parsed function
      71           8 :   _func_level_set = std::make_shared<SymFunction>();
      72          56 :   parsedFunctionSetup(_func_level_set,
      73           8 :                       _level_set,
      74             :                       "x,y,z",
      75             :                       getParam<std::vector<std::string>>("constant_names"),
      76             :                       getParam<std::vector<std::string>>("constant_expressions"),
      77             :                       comm());
      78           8 :   _func_params.resize(3);
      79             : 
      80             :   // Get the projection direction and normalize it
      81          24 :   auto proj_d = getParam<Point>("direction");
      82           8 :   if (MooseUtils::absoluteFuzzyEqual(proj_d.norm_sq(), 0))
      83           0 :     paramError("direction", "Should not be 0");
      84           8 :   _proj_dir = proj_d.unit();
      85           8 : }
      86             : 
      87             : std::unique_ptr<MeshBase>
      88           8 : ProjectSideSetOntoLevelSetGenerator::generate()
      89             : {
      90           8 :   auto mesh_uptr = std::move(_input);
      91           8 :   auto mesh = mesh_uptr.get();
      92           8 :   auto new_mesh = buildMeshBaseObject(3);
      93             : 
      94           8 :   if (!mesh_uptr->is_serial())
      95           0 :     paramError("input", "Input mesh must not be distributed");
      96             : 
      97             :   // Use a new subdomain ID to avoid an overlap if merging the two later
      98           8 :   const auto projection_block_id = MooseMeshUtils::getNextFreeSubdomainID(*mesh);
      99             : 
     100           8 :   BoundaryInfo & boundary_info = mesh->get_boundary_info();
     101           8 :   const auto bdry_side_list = boundary_info.build_side_list();
     102          16 :   const auto side_id = boundary_info.get_id_by_name(getParam<BoundaryName>("sideset"));
     103             : 
     104             :   // Loop over the sides on a given set
     105         424 :   for (const auto & [elem_id, side_i, bc_id] : bdry_side_list)
     106             :   {
     107         416 :     if (bc_id != side_id)
     108         208 :       continue;
     109             : 
     110             :     // Find the intersection with the surface (projection) for each node of the side
     111         208 :     std::unique_ptr<Elem> side_elem = mesh->elem_ptr(elem_id)->build_side_ptr(side_i);
     112         208 :     const auto n_nodes = side_elem->n_nodes();
     113         208 :     const auto side_type = side_elem->type();
     114         208 :     std::vector<Node *> new_nodes(n_nodes);
     115             : 
     116        1040 :     for (const auto i : make_range(n_nodes))
     117             :     {
     118         832 :       const auto start = side_elem->point(i);
     119         832 :       const auto end = start + _max_search_distance * _proj_dir;
     120             : 
     121         832 :       new_nodes[i] = new_mesh->add_point(pointPairLevelSetInterception(start, end));
     122             :     }
     123             : 
     124             :     // Build the same element as the side, on the surface, using the projected nodes
     125         208 :     if (side_type == libMesh::C0POLYGON)
     126           0 :       mooseError("Projection of polygonal sides is not supported at this time");
     127         208 :     std::unique_ptr<Elem> new_elem = Elem::build(side_type);
     128        1040 :     for (const auto i : make_range(n_nodes))
     129         832 :       new_elem->set_node(i, new_nodes[i]);
     130             : 
     131             :     // User-selected block name: same for all element types for now
     132         208 :     new_elem->subdomain_id() = projection_block_id;
     133             : 
     134         208 :     new_mesh->add_elem(std::move(new_elem));
     135         208 :   }
     136             : 
     137          16 :   new_mesh->subdomain_name(projection_block_id) = getParam<SubdomainName>("subdomain_name");
     138           8 :   new_mesh->unset_is_prepared();
     139          16 :   return dynamic_pointer_cast<MeshBase>(new_mesh);
     140           8 : }
     141             : 
     142             : Point
     143         832 : ProjectSideSetOntoLevelSetGenerator::pointPairLevelSetInterception(const Point & point1,
     144             :                                                                    const Point & point2)
     145             : {
     146         832 :   Real dist1 = levelSetEvaluator(point1);
     147         832 :   Real dist2 = levelSetEvaluator(point2);
     148             : 
     149         832 :   if (MooseUtils::absoluteFuzzyEqual(dist1, 0.0) || MooseUtils::absoluteFuzzyEqual(dist2, 0.0))
     150           0 :     mooseError("At least one of the two points are on the plane.");
     151         832 :   if (MooseUtils::absoluteFuzzyGreaterThan(dist1 * dist2, 0.0))
     152           0 :     mooseError("The two points are on the same side of the plane.");
     153             : 
     154         832 :   Point p1(point1);
     155         832 :   Point p2(point2);
     156         832 :   Real dist = abs(dist1) + abs(dist2);
     157         832 :   Point mid_point;
     158             : 
     159             :   // Bisection method to find midpoint
     160         832 :   unsigned int num_its = 0;
     161       51584 :   while (MooseUtils::absoluteFuzzyGreaterThan(dist, 0.0) && num_its < 1e3)
     162             :   {
     163       50752 :     num_its++;
     164       50752 :     mid_point = 0.5 * (p1 + p2);
     165       50752 :     const Real dist_mid = levelSetEvaluator(mid_point);
     166             :     // We do not need Fuzzy here as it will be covered by the while loop
     167       50752 :     if (dist_mid == 0.0)
     168           0 :       dist = 0.0;
     169       50752 :     else if (dist_mid * dist1 < 0.0)
     170             :     {
     171       32480 :       p2 = mid_point;
     172       32480 :       dist2 = levelSetEvaluator(p2);
     173       32480 :       dist = abs(dist1) + abs(dist2);
     174             :     }
     175             :     else
     176             :     {
     177       18272 :       p1 = mid_point;
     178       18272 :       dist1 = levelSetEvaluator(p1);
     179       18272 :       dist = abs(dist1) + abs(dist2);
     180             :     }
     181             :   }
     182             : 
     183         832 :   if (num_its == 1e3)
     184           0 :     mooseError("Projection failed for point " + Moose::stringify(point1) +
     185             :                ". Is the level set 'in front of' the sideset when following the projection "
     186             :                "direction?");
     187        1664 :   return mid_point;
     188             : }
     189             : 
     190             : Real
     191      103168 : ProjectSideSetOntoLevelSetGenerator::levelSetEvaluator(const Point & point)
     192             : {
     193      103168 :   _func_params[0] = point(0);
     194      103168 :   _func_params[1] = point(1);
     195      103168 :   _func_params[2] = point(2);
     196      309504 :   return evaluate(_func_level_set);
     197             : }