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 "BreakMeshByElementGenerator.h"
      11             : #include "CastUniquePointer.h"
      12             : #include "MooseMeshUtils.h"
      13             : 
      14             : #include "libmesh/partitioner.h"
      15             : 
      16             : registerMooseObject("MooseApp", BreakMeshByElementGenerator);
      17             : registerMooseObjectRenamed("MooseApp",
      18             :                            ExplodeMeshGenerator,
      19             :                            "05/18/2024 24:00",
      20             :                            BreakMeshByElementGenerator);
      21             : 
      22             : InputParameters
      23       28570 : BreakMeshByElementGenerator::validParams()
      24             : {
      25       28570 :   InputParameters params = MeshGenerator::validParams();
      26       28570 :   params.addClassDescription("Break all element-element interfaces in the specified subdomains.");
      27       28570 :   params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
      28       28570 :   params.addParam<std::vector<SubdomainID>>("subdomains", "The list of subdomain IDs to explode.");
      29       28570 :   params.addRequiredParam<BoundaryName>(
      30             :       "interface_name", "The boundary name containing all broken element-element interfaces.");
      31       28570 :   return params;
      32           0 : }
      33             : 
      34          20 : BreakMeshByElementGenerator::BreakMeshByElementGenerator(const InputParameters & parameters)
      35             :   : MeshGenerator(parameters),
      36          20 :     _input(getMesh("input")),
      37          20 :     _subdomains(getParam<std::vector<SubdomainID>>("subdomains")),
      38          40 :     _interface_name(getParam<BoundaryName>("interface_name"))
      39             : {
      40          20 : }
      41             : 
      42             : std::unique_ptr<MeshBase>
      43          20 : BreakMeshByElementGenerator::generate()
      44             : {
      45          20 :   std::unique_ptr<MeshBase> mesh = std::move(_input);
      46             : 
      47             :   // check that the subdomain IDs exist in the mesh
      48          52 :   for (const auto & id : _subdomains)
      49          36 :     if (!MooseMeshUtils::hasSubdomainID(*mesh, id))
      50           4 :       paramError("subdomains", "The block ID '", id, "' was not found in the mesh");
      51             : 
      52          16 :   BoundaryInfo & boundary_info = mesh->get_boundary_info();
      53          16 :   if (boundary_info.get_id_by_name(_interface_name) != Moose::INVALID_BOUNDARY_ID)
      54           0 :     paramError("interface_name", "The specified interface name already exits in the mesh.");
      55             : 
      56          16 :   const auto node_to_elem_map = buildSubdomainRestrictedNodeToElemMap(mesh, _subdomains);
      57             : 
      58          16 :   duplicateNodes(mesh, node_to_elem_map);
      59             : 
      60          16 :   createInterface(*mesh, node_to_elem_map);
      61             : 
      62          16 :   Partitioner::set_node_processor_ids(*mesh);
      63             : 
      64          32 :   return dynamic_pointer_cast<MeshBase>(mesh);
      65          16 : }
      66             : 
      67             : BreakMeshByElementGenerator::NodeToElemMapType
      68          16 : BreakMeshByElementGenerator::buildSubdomainRestrictedNodeToElemMap(
      69             :     std::unique_ptr<MeshBase> & mesh, const std::vector<SubdomainID> & subdomains) const
      70             : {
      71          16 :   NodeToElemMapType node_to_elem_map;
      72         784 :   for (const auto & elem : mesh->active_element_ptr_range())
      73             :   {
      74             :     // Skip if the element is not in the specified subdomains
      75         384 :     if (std::find(subdomains.begin(), subdomains.end(), elem->subdomain_id()) == subdomains.end())
      76         192 :       continue;
      77             : 
      78         192 :     std::set<const Elem *> neighbors;
      79         192 :     elem->find_point_neighbors(neighbors);
      80             : 
      81        1472 :     for (auto n : make_range(elem->n_nodes()))
      82             :     {
      83             :       // if ANY neighboring element that contains this node is not in the specified subdomains,
      84             :       // don't add this node to the map, i.e. don't split this node.
      85        1280 :       bool should_duplicate = true;
      86       12253 :       for (auto neighbor : neighbors)
      87       15616 :         if (neighbor->contains_point(elem->node_ref(n)) &&
      88        4323 :             std::find(subdomains.begin(), subdomains.end(), neighbor->subdomain_id()) ==
      89       15616 :                 subdomains.end())
      90             :         {
      91         320 :           should_duplicate = false;
      92         320 :           break;
      93             :         }
      94             : 
      95        1280 :       if (should_duplicate)
      96         960 :         node_to_elem_map[elem->node_id(n)].insert(elem->id());
      97             :     }
      98         208 :   }
      99             : 
     100          16 :   return node_to_elem_map;
     101           0 : }
     102             : 
     103             : void
     104          16 : BreakMeshByElementGenerator::duplicateNodes(std::unique_ptr<MeshBase> & mesh,
     105             :                                             const NodeToElemMapType & node_to_elem_map) const
     106             : {
     107         336 :   for (const auto & [node_id, connected_elem_ids] : node_to_elem_map)
     108        1280 :     for (auto & connected_elem_id : connected_elem_ids)
     109         960 :       if (connected_elem_id != *connected_elem_ids.begin())
     110         640 :         duplicateNode(mesh, mesh->elem_ptr(connected_elem_id), mesh->node_ptr(node_id));
     111          16 : }
     112             : 
     113             : void
     114         640 : BreakMeshByElementGenerator::duplicateNode(std::unique_ptr<MeshBase> & mesh,
     115             :                                            Elem * elem,
     116             :                                            const Node * node) const
     117             : {
     118         640 :   std::unique_ptr<Node> new_node = Node::build(*node, Node::invalid_id);
     119         640 :   new_node->processor_id() = elem->processor_id();
     120         640 :   Node * added_node = mesh->add_node(std::move(new_node));
     121        2400 :   for (const auto j : elem->node_index_range())
     122        2400 :     if (elem->node_id(j) == node->id())
     123             :     {
     124         640 :       elem->set_node(j, added_node);
     125         640 :       break;
     126             :     }
     127             : 
     128             :   // Add boundary info to the new node
     129         640 :   BoundaryInfo & boundary_info = mesh->get_boundary_info();
     130         640 :   std::vector<boundary_id_type> node_boundary_ids;
     131         640 :   boundary_info.boundary_ids(node, node_boundary_ids);
     132         640 :   boundary_info.add_node(added_node, node_boundary_ids);
     133         640 : }
     134             : 
     135             : void
     136          16 : BreakMeshByElementGenerator::createInterface(MeshBase & mesh,
     137             :                                              const NodeToElemMapType & node_to_elem_map) const
     138             : {
     139          16 :   BoundaryInfo & boundary_info = mesh.get_boundary_info();
     140          16 :   const auto & existing_boundary_ids = boundary_info.get_boundary_ids();
     141             :   boundary_id_type interface_id =
     142          16 :       existing_boundary_ids.empty() ? 0 : *existing_boundary_ids.rbegin() + 1;
     143          16 :   boundary_info.sideset_name(interface_id) = _interface_name;
     144             : 
     145          16 :   std::set<std::pair<dof_id_type, unsigned int>> sides_to_add;
     146             : 
     147         336 :   for (const auto & node_to_elems : node_to_elem_map)
     148        1280 :     for (const auto & elem_id_i : node_to_elems.second)
     149             :     {
     150         960 :       Elem * elem_i = mesh.elem_ptr(elem_id_i);
     151        4880 :       for (const auto & elem_id_j : node_to_elems.second)
     152             :       {
     153        3920 :         Elem * elem_j = mesh.elem_ptr(elem_id_j);
     154        3920 :         if (elem_i != elem_j && elem_id_i < elem_id_j && elem_i->has_neighbor(elem_j))
     155         872 :           sides_to_add.insert(std::make_pair(elem_id_i, elem_i->which_neighbor_am_i(elem_j)));
     156             :       }
     157             :     }
     158             : 
     159         320 :   for (const auto & [elem_id, side] : sides_to_add)
     160         304 :     boundary_info.add_side(elem_id, side, interface_id);
     161          16 : }