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 "GapLineMeshGenerator.h"
      11             : 
      12             : #include "MooseMeshUtils.h"
      13             : #include "MooseUtils.h"
      14             : #include "GeometryUtils.h"
      15             : 
      16             : #include "libmesh/mesh_serializer.h"
      17             : #include "libmesh/mesh_triangle_holes.h"
      18             : #include "libmesh/poly2tri_triangulator.h"
      19             : 
      20             : registerMooseObject("ReactorApp", GapLineMeshGenerator);
      21             : 
      22             : InputParameters
      23          60 : GapLineMeshGenerator::validParams()
      24             : {
      25          60 :   InputParameters params = PolygonMeshGeneratorBase::validParams();
      26             : 
      27         120 :   params.addRequiredParam<MeshGeneratorName>("input", "The input mesh to create the gap based on.");
      28             : 
      29         120 :   params.addRequiredParam<Real>("thickness", "The thickness of the gap to be created.");
      30             : 
      31         120 :   MooseEnum gap_direction("OUTWARD INWARD", "OUTWARD");
      32             : 
      33         120 :   params.addParam<MooseEnum>("gap_direction",
      34             :                              gap_direction,
      35             :                              "In which direction the gap is created with respect to the side "
      36             :                              "normal of the elements along the boundary of the input mesh.");
      37             : 
      38          60 :   params.addParam<std::vector<boundary_id_type>>(
      39             :       "boundary_ids",
      40          60 :       std::vector<boundary_id_type>(),
      41             :       "The boundary IDs around which the gap will be created.");
      42             : 
      43         120 :   params.addParam<Real>("max_elem_size", "The maximum element size for the generated gap mesh.");
      44             : 
      45          60 :   params.addClassDescription(
      46             :       "Generates a polyline mesh that is based on an input 2D-XY mesh. The 2D-XY mesh needs to be "
      47             :       "a "
      48             :       "connected mesh with only one outer boundary manifold. The polyline mesh generated along "
      49             :       "with the boundary of the input mesh form an unmeshed gap with a specified thickness.");
      50             : 
      51          60 :   return params;
      52          60 : }
      53             : 
      54          30 : GapLineMeshGenerator::GapLineMeshGenerator(const InputParameters & parameters)
      55             :   : PolygonMeshGeneratorBase(parameters),
      56          30 :     _input(getMesh("input")),
      57          60 :     _thickness(getParam<Real>("thickness")),
      58          60 :     _gap_direction(getParam<MooseEnum>("gap_direction").template getEnum<GapDirection>()),
      59          90 :     _boundary_ids(getParam<std::vector<boundary_id_type>>("boundary_ids"))
      60             : {
      61          30 : }
      62             : 
      63             : std::unique_ptr<MeshBase>
      64          30 : GapLineMeshGenerator::generate()
      65             : {
      66             :   // Put the boundary mesh in a local pointer
      67             :   std::unique_ptr<UnstructuredMesh> mesh =
      68          30 :       dynamic_pointer_cast<UnstructuredMesh>(std::move(_input));
      69             : 
      70          30 :   std::set<std::size_t> mesh_bdry_ids(_boundary_ids.begin(), _boundary_ids.end());
      71             :   // MeshedHole is a good tool to extract and sort boundary points
      72          60 :   TriangulatorInterface::MeshedHole bdry_mh(*mesh, mesh_bdry_ids);
      73             : 
      74             :   // Reduce the point list to only contain vertices
      75             :   std::vector<Point> reduced_pts_list;
      76         390 :   for (const auto i : make_range(bdry_mh.n_points()))
      77             :   {
      78         360 :     if (!geom_utils::arePointsColinear(
      79         360 :             bdry_mh.point((i - 1 + bdry_mh.n_points()) % bdry_mh.n_points()),
      80         360 :             bdry_mh.point(i),
      81         720 :             bdry_mh.point((i + 1) % bdry_mh.n_points())))
      82         720 :       reduced_pts_list.push_back(bdry_mh.point(i));
      83             :   }
      84             :   // Here we need a method to generate the outward normals of each external side
      85          30 :   auto ply_mesh = buildMeshBaseObject();
      86             : 
      87          30 :   MooseMeshUtils::buildPolyLineMesh(*ply_mesh,
      88             :                                     reduced_pts_list,
      89             :                                     /*loop*/ true,
      90          30 :                                     BoundaryName(),
      91          30 :                                     BoundaryName(),
      92          60 :                                     std::vector<unsigned int>({1}));
      93             : 
      94             :   std::unique_ptr<UnstructuredMesh> ply_mesh_u =
      95          30 :       dynamic_pointer_cast<UnstructuredMesh>(std::move(ply_mesh));
      96             : 
      97             :   // Generate a very simple triangulation mesh so that we can get the outward normal vectors
      98          30 :   libMesh::Poly2TriTriangulator poly2tri(*ply_mesh_u);
      99          30 :   poly2tri.triangulation_type() = libMesh::TriangulatorInterface::PSLG;
     100             : 
     101          30 :   poly2tri.set_interpolate_boundary_points(0);
     102             :   poly2tri.set_refine_boundary_allowed(false);
     103             :   poly2tri.set_verify_hole_boundaries(false);
     104          30 :   poly2tri.desired_area() = 0;
     105          30 :   poly2tri.minimum_angle() = 0; // Not yet supported
     106          30 :   poly2tri.smooth_after_generating() = false;
     107          30 :   poly2tri.triangulate();
     108             : 
     109             :   // We need to serialize the mesh for next steps
     110          30 :   libMesh::MeshSerializer serial(*ply_mesh_u);
     111             :   // The mesh now only contains one side set that corresponds to the outer boundary with an ID of 0
     112          30 :   auto bdry_list(ply_mesh_u->get_boundary_info().build_side_list());
     113             : 
     114             :   // For each vertex, the shifting direction to form the gap is defined by the normal vectors of the
     115             :   // two sides that contain the vertex
     116             :   // We gather the normals for each node on the boundary here, which are all vertices because of the
     117             :   // pre-selection of nodes
     118             :   std::map<dof_id_type, std::vector<Point>> node_normal_map;
     119         390 :   for (const auto & bside : bdry_list)
     120             :   {
     121         360 :     const auto & side = ply_mesh_u->elem_ptr(std::get<0>(bside))->side_ptr(std::get<1>(bside));
     122             :     const Point side_normal =
     123         360 :         ply_mesh_u->elem_ptr(std::get<0>(bside))->side_vertex_average_normal(std::get<1>(bside));
     124         360 :     if (node_normal_map.count(side->node_ptr(0)->id()))
     125         218 :       node_normal_map[side->node_ptr(0)->id()].push_back(side_normal);
     126             :     else
     127         142 :       node_normal_map[side->node_ptr(0)->id()] = {side_normal};
     128         360 :     if (node_normal_map.count(side->node_ptr(1)->id()))
     129         142 :       node_normal_map[side->node_ptr(1)->id()].push_back(side_normal);
     130             :     else
     131         218 :       node_normal_map[side->node_ptr(1)->id()] = {side_normal};
     132         360 :   }
     133             : 
     134          30 :   std::vector<Point> mod_reduced_pts_list(reduced_pts_list);
     135         390 :   for (const auto & [node_id, normal_vecs] : node_normal_map)
     136             :   {
     137             :     mooseAssert(normal_vecs.size() == 2,
     138             :                 "Each vertex should be connected to exactly two sides in a polygon.");
     139             : 
     140         360 :     const Point original_pt = *(ply_mesh_u->node_ptr(node_id));
     141             :     // Form an average normal at the vertex from the two connected sides' normals
     142         360 :     const Point move_dir = (normal_vecs.front() + normal_vecs.back()).unit() *
     143         360 :                            ((_gap_direction == GapDirection::OUTWARD) ? 1.0 : -1.0);
     144             :     // Consider four points of interest to determine the moving distance
     145             :     // 1. the vertex point
     146             :     // 2. point along normal_vecs.front() from the vertex with a distance _thickness
     147             :     // 3. point along normal_vecs.back() from the vertex with a distance _thickness
     148             :     // 4. point along move_dir from the vertex with a distance mov_dist
     149             :     // The four points form a kite shape with its symmetry axis along 1-4 direction
     150             :     // Angle 1-2-4 and angle 1-3-4 are right angles
     151             :     // Angle 2-1-3 (i.e., theta) can be calculated using the interior product of
     152             :     // v1 = normal_vecs.front() and v2 = normal_vecs.back():
     153             :     // cos(theta) = (v1 . v2) / (|v1| * |v2|)
     154             :     // Because of the right angles, the distance between 1 and 4 can be calculated as:
     155             :     // mov_dist = _thickness / cos(theta/2)
     156             :     // and cos(theta/2) = sqrt((1 + cos(theta)) / 2)
     157             :     const Real mov_dist =
     158         360 :         _thickness /
     159         360 :         std::sqrt((1.0 + (normal_vecs.front() * normal_vecs.back()) /
     160         360 :                              (normal_vecs.front().norm() * normal_vecs.back().norm())) /
     161             :                   2.0);
     162             : 
     163             :     mooseAssert(std::count(reduced_pts_list.begin(), reduced_pts_list.end(), original_pt) == 1,
     164             :                 "The original point should be found exactly once in the reduced points list.");
     165         360 :     mod_reduced_pts_list[std::distance(
     166             :         reduced_pts_list.begin(),
     167         360 :         std::find(reduced_pts_list.begin(), reduced_pts_list.end(), original_pt))] =
     168             :         original_pt + move_dir * mov_dist;
     169             :   }
     170             : 
     171             :   // To ensure no overlapping, we need to check set of four points
     172             :   // p1 and p2 should be the pair of points before shifting the side
     173             :   // p3 and p4 should be the pair of points after shifting a side
     174         375 :   for (const auto & i_node_1 : make_range(mod_reduced_pts_list.size()))
     175             :   {
     176             :     const Point & p1 = reduced_pts_list[i_node_1];
     177             :     const Point & p2 = mod_reduced_pts_list[i_node_1];
     178        4504 :     for (const auto & i_node_2 : make_range(mod_reduced_pts_list.size()))
     179             :     {
     180        4159 :       if (i_node_2 == i_node_1 || (i_node_2 + 1) % mod_reduced_pts_list.size() == i_node_1)
     181             :         continue;
     182             :       const Point & p3 = mod_reduced_pts_list[i_node_2];
     183             :       const Point & p4 = mod_reduced_pts_list[(i_node_2 + 1) % mod_reduced_pts_list.size()];
     184        3465 :       if (geom_utils::segmentsIntersect(p1, p2, p3, p4))
     185           2 :         paramError("thickness",
     186             :                    "The thickness is so large that the mesh is tangled because the offset nodes "
     187             :                    "are no longer in the same order when following the original boundary. Please "
     188             :                    "reduce the thickness value.");
     189             :     }
     190             :   }
     191             : 
     192          28 :   auto ply_mesh_2 = buildMeshBaseObject();
     193             : 
     194          56 :   if (isParamValid("max_elem_size"))
     195           7 :     MooseMeshUtils::buildPolyLineMesh(*ply_mesh_2,
     196             :                                       mod_reduced_pts_list,
     197             :                                       true,
     198           7 :                                       BoundaryName(),
     199           7 :                                       BoundaryName(),
     200          14 :                                       getParam<Real>("max_elem_size"));
     201             :   else
     202          21 :     MooseMeshUtils::buildPolyLineMesh(*ply_mesh_2,
     203             :                                       mod_reduced_pts_list,
     204             :                                       true,
     205          21 :                                       BoundaryName(),
     206          21 :                                       BoundaryName(),
     207          42 :                                       std::vector<unsigned int>({1}));
     208             : 
     209          28 :   return ply_mesh_2;
     210          84 : }