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 "PlaneIDMeshGenerator.h"
      11             : #include "CastUniquePointer.h"
      12             : 
      13             : #include "libmesh/elem.h"
      14             : 
      15             : registerMooseObject("MooseApp", PlaneIDMeshGenerator);
      16             : 
      17             : InputParameters
      18       14381 : PlaneIDMeshGenerator::validParams()
      19             : {
      20       14381 :   InputParameters params = MeshGenerator::validParams();
      21             : 
      22       14381 :   params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
      23       14381 :   params.addRequiredParam<std::vector<Real>>(
      24             :       "plane_coordinates",
      25             :       "Coordinates of planes along the axis. The origin are at x/y/z=0 depending on the axis");
      26       14381 :   params.addParam<std::vector<unsigned int>>("num_ids_per_plane", "Number of unique ids per plane");
      27       14381 :   MooseEnum plane_axis("x y z", "z");
      28       14381 :   params.addParam<MooseEnum>("plane_axis", plane_axis, "Axis of plane");
      29       14381 :   params.addRequiredParam<std::string>("id_name", "Name of extra integer ID set");
      30       14381 :   params.addParam<Real>("tolerance", 1.0E-4, "Tolerance for plane coordinate check");
      31       14381 :   params.addClassDescription("Adds an extra element integer that identifies planes in a mesh.");
      32       28762 :   return params;
      33       14381 : }
      34             : 
      35          64 : PlaneIDMeshGenerator::PlaneIDMeshGenerator(const InputParameters & parameters)
      36             :   : MeshGenerator(parameters),
      37          64 :     _input(getMesh("input")),
      38          64 :     _axis_index(getParam<MooseEnum>("plane_axis")),
      39         192 :     _element_id_name(getParam<std::string>("id_name"))
      40             : {
      41          64 :   if (!isParamValid("num_ids_per_plane"))
      42           0 :     _planes = getParam<std::vector<Real>>("plane_coordinates");
      43             :   else
      44             :   {
      45          64 :     _planes.clear();
      46          64 :     std::vector<Real> base_planes = getParam<std::vector<Real>>("plane_coordinates");
      47          64 :     std::vector<unsigned int> sublayers = getParam<std::vector<unsigned int>>("num_ids_per_plane");
      48          64 :     if (base_planes.size() != sublayers.size() + 1)
      49           4 :       paramError("plane_coordinates",
      50             :                  "Sizes of 'plane_coordinates' and 'num_ids_per_plane' disagree");
      51          60 :     _planes.push_back(base_planes[0]);
      52         184 :     for (unsigned int i = 0; i < sublayers.size(); ++i)
      53             :     {
      54         128 :       Real layer_size = (base_planes[i + 1] - base_planes[i]) / (Real)sublayers[i];
      55         128 :       if (MooseUtils::absoluteFuzzyGreaterEqual(base_planes[i], base_planes[i + 1]))
      56           4 :         paramError("plane_coordinates", "Plane coordinates must be in increasing order");
      57             : 
      58         300 :       for (unsigned int j = 0; j < sublayers[i]; ++j)
      59         176 :         _planes.push_back(_planes.back() + layer_size);
      60             :     }
      61          56 :   }
      62          56 :   if (_planes.size() < 2)
      63           4 :     paramError("plane_coordinates", "Size of 'plane_coordinates' should be at least two");
      64          52 : }
      65             : 
      66             : std::unique_ptr<MeshBase>
      67          52 : PlaneIDMeshGenerator::generate()
      68             : {
      69          52 :   std::unique_ptr<MeshBase> mesh = std::move(_input);
      70          52 :   if (mesh->mesh_dimension() < _axis_index + 1)
      71           4 :     paramError("plane_axis", "Plane axis must be contained in the mesh. Check mesh dimension");
      72          48 :   unsigned int extra_id_index = 0;
      73          48 :   if (!mesh->has_elem_integer(_element_id_name))
      74          48 :     extra_id_index = mesh->add_elem_integer(_element_id_name);
      75             :   else
      76             :   {
      77           0 :     extra_id_index = mesh->get_elem_integer_index(_element_id_name);
      78           0 :     paramWarning(
      79           0 :         "id_name", "An element integer with the name '", _element_id_name, "' already exists");
      80             :   }
      81             : 
      82             :   // Check that the planes do not cut the elements beyond the acceptable tolerance
      83          48 :   const Real tol = getParam<Real>("tolerance");
      84        7574 :   for (auto & elem : mesh->active_element_ptr_range())
      85             :   {
      86        3771 :     const int layer_id = getPlaneID(elem->vertex_average());
      87       31299 :     for (unsigned int i = 0; i < elem->n_nodes(); ++i)
      88             :     {
      89       27536 :       const Point & p = elem->point(i) - (elem->point(i) - elem->vertex_average()) * tol;
      90       27536 :       if (getPlaneID(p) != layer_id)
      91           4 :         mooseError("Element at ", elem->vertex_average(), " is cut by the planes");
      92             :     }
      93        3763 :     elem->set_extra_integer(extra_id_index, layer_id);
      94          40 :   }
      95          80 :   return dynamic_pointer_cast<MeshBase>(mesh);
      96          40 : }
      97             : 
      98             : int
      99       31307 : PlaneIDMeshGenerator::getPlaneID(const Point & p) const
     100             : {
     101             :   // check the input point is located within the plane range
     102       31307 :   if (p(_axis_index) < _planes[0] || p(_axis_index) > _planes.back())
     103           4 :     mooseError("The planes do not cover element at ", p);
     104       31303 :   int id = 0;
     105             :   // looping over each plane to find a plane ID corresponding to the input point
     106       78157 :   for (unsigned int layer_id = 0; layer_id < _planes.size() - 1; ++layer_id)
     107             :   {
     108             :     // check the input point is located between _planes[layer_id] and _planes[layer_id + 1]
     109       78157 :     if (_planes[layer_id] < p(_axis_index) && _planes[layer_id + 1] >= p(_axis_index))
     110             :     {
     111       31303 :       id = layer_id;
     112       31303 :       break;
     113             :     }
     114             :   }
     115       31303 :   return id;
     116             : }