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 "LowerDBlockFromSidesetGenerator.h"
      11             : #include "InputParameters.h"
      12             : #include "MooseTypes.h"
      13             : #include "CastUniquePointer.h"
      14             : #include "MooseMeshUtils.h"
      15             : 
      16             : #include "libmesh/distributed_mesh.h"
      17             : #include "libmesh/elem.h"
      18             : #include "libmesh/parallel_elem.h"
      19             : #include "libmesh/parallel_node.h"
      20             : #include "libmesh/compare_elems_by_level.h"
      21             : #include "libmesh/mesh_communication.h"
      22             : 
      23             : #include "timpi/parallel_sync.h"
      24             : 
      25             : #include <set>
      26             : #include <typeinfo>
      27             : 
      28             : registerMooseObject("MooseApp", LowerDBlockFromSidesetGenerator);
      29             : 
      30             : InputParameters
      31       18503 : LowerDBlockFromSidesetGenerator::validParams()
      32             : {
      33       18503 :   InputParameters params = MeshGenerator::validParams();
      34             : 
      35       18503 :   params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
      36       18503 :   params.addParam<SubdomainID>("new_block_id", "The lower dimensional block id to create");
      37       18503 :   params.addParam<SubdomainName>("new_block_name",
      38             :                                  "The lower dimensional block name to create (optional)");
      39       18503 :   params.addRequiredParam<std::vector<BoundaryName>>(
      40             :       "sidesets", "The sidesets from which to create the new block");
      41             : 
      42       18503 :   params.addClassDescription("Adds lower dimensional elements on the specified sidesets.");
      43             : 
      44       18503 :   return params;
      45           0 : }
      46             : 
      47        2119 : LowerDBlockFromSidesetGenerator::LowerDBlockFromSidesetGenerator(const InputParameters & parameters)
      48             :   : MeshGenerator(parameters),
      49        2119 :     _input(getMesh("input")),
      50        4238 :     _sideset_names(getParam<std::vector<BoundaryName>>("sidesets"))
      51             : {
      52        2119 : }
      53             : 
      54             : // Used to temporarily store information about which lower-dimensional
      55             : // sides to add and what subdomain id to use for the added sides.
      56             : struct ElemSideDouble
      57             : {
      58       19110 :   ElemSideDouble(Elem * elem_in, unsigned short int side_in) : elem(elem_in), side(side_in) {}
      59             : 
      60             :   Elem * elem;
      61             :   unsigned short int side;
      62             : };
      63             : 
      64             : std::unique_ptr<MeshBase>
      65        1986 : LowerDBlockFromSidesetGenerator::generate()
      66             : {
      67        1986 :   std::unique_ptr<MeshBase> mesh = std::move(_input);
      68             : 
      69             :   // Generate a new block id if one isn't supplied.
      70        3972 :   SubdomainID new_block_id = isParamValid("new_block_id")
      71        3972 :                                  ? getParam<SubdomainID>("new_block_id")
      72        1986 :                                  : MooseMeshUtils::getNextFreeSubdomainID(*mesh);
      73             : 
      74             :   // Make sure our boundary info and parallel counts are setup
      75        1986 :   if (!mesh->is_prepared())
      76             :   {
      77         384 :     const bool allow_remote_element_removal = mesh->allow_remote_element_removal();
      78             :     // We want all of our boundary elements available, so avoid removing them if they haven't
      79             :     // already been so
      80         384 :     mesh->allow_remote_element_removal(false);
      81         384 :     mesh->prepare_for_use();
      82         384 :     mesh->allow_remote_element_removal(allow_remote_element_removal);
      83             :   }
      84             : 
      85             :   // Check that the sidesets are present in the mesh
      86        4132 :   for (const auto & sideset : _sideset_names)
      87        2150 :     if (!MooseMeshUtils::hasBoundaryName(*mesh, sideset))
      88           4 :       paramError("sidesets", "The sideset '", sideset, "' was not found within the mesh");
      89             : 
      90        1982 :   auto sideset_ids = MooseMeshUtils::getBoundaryIDs(*mesh, _sideset_names, true);
      91        1982 :   std::set<boundary_id_type> sidesets(sideset_ids.begin(), sideset_ids.end());
      92        1982 :   auto side_list = mesh->get_boundary_info().build_side_list();
      93        1982 :   if (!mesh->is_serial() && mesh->comm().size() > 1)
      94             :   {
      95          10 :     std::vector<Elem *> elements_to_send;
      96          10 :     unsigned short i_need_boundary_elems = 0;
      97         172 :     for (const auto & [elem_id, side, bc_id] : side_list)
      98             :     {
      99         162 :       libmesh_ignore(side);
     100         162 :       if (sidesets.count(bc_id))
     101             :       {
     102             :         // Whether we have this boundary information through our locally owned element or a ghosted
     103             :         // element, we'll need the boundary elements for parallel consistent addition
     104          92 :         i_need_boundary_elems = 1;
     105          92 :         auto * elem = mesh->elem_ptr(elem_id);
     106          92 :         if (elem->processor_id() == mesh->processor_id())
     107          80 :           elements_to_send.push_back(elem);
     108             :       }
     109             :     }
     110             : 
     111             :     std::set<const Elem *, libMesh::CompareElemIdsByLevel> connected_elements(
     112          10 :         elements_to_send.begin(), elements_to_send.end());
     113          10 :     std::set<const Node *> connected_nodes;
     114          10 :     reconnect_nodes(connected_elements, connected_nodes);
     115          10 :     std::set<dof_id_type> connected_node_ids;
     116         178 :     for (auto * nd : connected_nodes)
     117         168 :       connected_node_ids.insert(nd->id());
     118             : 
     119          10 :     std::vector<unsigned short> need_boundary_elems(mesh->comm().size());
     120          10 :     mesh->comm().allgather(i_need_boundary_elems, need_boundary_elems);
     121          10 :     std::unordered_map<processor_id_type, decltype(elements_to_send)> push_element_data;
     122          10 :     std::unordered_map<processor_id_type, decltype(connected_nodes)> push_node_data;
     123             : 
     124          30 :     for (const auto pid : index_range(mesh->comm()))
     125             :       // Don't need to send to self
     126          20 :       if (pid != mesh->processor_id() && need_boundary_elems[pid])
     127             :       {
     128          10 :         if (elements_to_send.size())
     129          10 :           push_element_data[pid] = elements_to_send;
     130          10 :         if (connected_nodes.size())
     131          10 :           push_node_data[pid] = connected_nodes;
     132             :       }
     133             : 
     134          10 :     auto node_action_functor = [](processor_id_type, const auto &)
     135             :     {
     136             :       // Node packing specialization already has unpacked node into mesh, so nothing to do
     137          10 :     };
     138          20 :     Parallel::push_parallel_packed_range(
     139          10 :         mesh->comm(), push_node_data, mesh.get(), node_action_functor);
     140          10 :     auto elem_action_functor = [](processor_id_type, const auto &)
     141             :     {
     142             :       // Elem packing specialization already has unpacked elem into mesh, so nothing to do
     143          10 :     };
     144          20 :     TIMPI::push_parallel_packed_range(
     145          10 :         mesh->comm(), push_element_data, mesh.get(), elem_action_functor);
     146             : 
     147             :     // now that we've gathered everything, we need to rebuild the side list
     148          10 :     side_list = mesh->get_boundary_info().build_side_list();
     149          10 :   }
     150             : 
     151        1982 :   std::vector<std::pair<dof_id_type, ElemSideDouble>> element_sides_on_boundary;
     152        1982 :   dof_id_type counter = 0;
     153      127604 :   for (const auto & triple : side_list)
     154      125626 :     if (sidesets.count(std::get<2>(triple)))
     155             :     {
     156       19114 :       if (auto elem = mesh->query_elem_ptr(std::get<0>(triple)))
     157             :       {
     158       19114 :         if (!elem->active())
     159           4 :           mooseError(
     160             :               "Only active, level 0 elements can be made interior parents of new level 0 lower-d "
     161             :               "elements. Make sure that ",
     162           4 :               type(),
     163             :               "s are run before any refinement generators");
     164       19110 :         element_sides_on_boundary.push_back(
     165       38220 :             std::make_pair(counter, ElemSideDouble(elem, std::get<1>(triple))));
     166             :       }
     167       19110 :       ++counter;
     168             :     }
     169             : 
     170        1978 :   dof_id_type max_elem_id = mesh->max_elem_id();
     171        1978 :   unique_id_type max_unique_id = mesh->parallel_max_unique_id();
     172             : 
     173             :   // Making an important assumption that at least our boundary elements are the same on all
     174             :   // processes even in distributed mesh mode (this is reliant on the correct ghosting functors
     175             :   // existing on the mesh)
     176       21088 :   for (auto & [i, elem_side] : element_sides_on_boundary)
     177             :   {
     178       19110 :     Elem * elem = elem_side.elem;
     179             : 
     180       19110 :     const auto side = elem_side.side;
     181             : 
     182             :     // Build a non-proxy element from this side.
     183       19110 :     std::unique_ptr<Elem> side_elem(elem->build_side_ptr(side));
     184             : 
     185             :     // The side will be added with the same processor id as the parent.
     186       19110 :     side_elem->processor_id() = elem->processor_id();
     187             : 
     188             :     // Add subdomain ID
     189       19110 :     side_elem->subdomain_id() = new_block_id;
     190             : 
     191             :     // Also assign the side's interior parent, so it is always
     192             :     // easy to figure out the Elem we came from.
     193       19110 :     side_elem->set_interior_parent(elem);
     194             : 
     195             :     // Add id
     196       19110 :     side_elem->set_id(max_elem_id + i);
     197       19110 :     side_elem->set_unique_id(max_unique_id + i);
     198             : 
     199             :     // Finally, add the lower-dimensional element to the Mesh.
     200       19110 :     mesh->add_elem(side_elem.release());
     201       19110 :   };
     202             : 
     203             :   // Assign block name, if provided
     204        1978 :   if (isParamValid("new_block_name"))
     205        1339 :     mesh->subdomain_name(new_block_id) = getParam<SubdomainName>("new_block_name");
     206             : 
     207        1978 :   const bool skip_partitioning_old = mesh->skip_partitioning();
     208        1978 :   mesh->skip_partitioning(true);
     209        1978 :   mesh->prepare_for_use();
     210        1978 :   mesh->skip_partitioning(skip_partitioning_old);
     211             : 
     212        3956 :   return mesh;
     213        1978 : }