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 6232 : BreakMeshByElementGenerator::validParams()
24 : {
25 6232 : InputParameters params = MeshGenerator::validParams();
26 12464 : params.addClassDescription("Break all element-element interfaces in the specified subdomains.");
27 24928 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
28 18696 : params.addParam<std::vector<SubdomainID>>(
29 : "subdomains",
30 12464 : std::vector<SubdomainID>(),
31 : "The list of subdomain IDs to explode. Leave unset to explode all subdomains.");
32 24928 : params.addParam<BoundaryName>(
33 : "interface_name",
34 : "element_boundaries",
35 : "The boundary name containing all broken element-element interfaces.");
36 24928 : params.addRangeCheckedParam<unsigned int>(
37 : "interface_sides",
38 12464 : 1,
39 : "interface_sides<3",
40 : "Whether to add no interface boundary, a 1-sided boundary (facing from lower to higher "
41 : "element id), or a 2-sided boundary");
42 6232 : return params;
43 0 : }
44 :
45 55 : BreakMeshByElementGenerator::BreakMeshByElementGenerator(const InputParameters & parameters)
46 : : MeshGenerator(parameters),
47 55 : _input(getMesh("input")),
48 110 : _subdomains(getParam<std::vector<SubdomainID>>("subdomains")),
49 110 : _interface_name(getParam<BoundaryName>("interface_name")),
50 165 : _interface_sides(getParam<unsigned int>("interface_sides"))
51 : {
52 55 : }
53 :
54 : std::unique_ptr<MeshBase>
55 52 : BreakMeshByElementGenerator::generate()
56 : {
57 52 : std::unique_ptr<MeshBase> mesh = std::move(_input);
58 :
59 52 : if (!mesh->is_prepared())
60 19 : mesh->prepare_for_use();
61 :
62 : // check that the subdomain IDs exist in the mesh
63 106 : for (const auto & id : _subdomains)
64 57 : if (!MooseMeshUtils::hasSubdomainID(*mesh, id))
65 6 : paramError("subdomains", "The block ID '", id, "' was not found in the mesh");
66 :
67 49 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
68 98 : if (_interface_sides &&
69 49 : boundary_info.get_id_by_name(_interface_name) != Moose::INVALID_BOUNDARY_ID)
70 0 : paramError("interface_name", "The specified interface name already exists in the mesh.");
71 :
72 49 : const auto node_to_elem_map = buildSubdomainRestrictedNodeToElemMap(mesh, _subdomains);
73 :
74 49 : duplicateNodes(mesh, node_to_elem_map);
75 :
76 49 : createInterface(*mesh, node_to_elem_map);
77 :
78 49 : Partitioner::set_node_processor_ids(*mesh);
79 :
80 : // New nodes were added with non-contiguous IDs (DistributedMesh's
81 : // unpartitioned counter strides by n_procs+1). Renumber so that
82 : // max_node_id == n_nodes, which is required by non-parallel output formats.
83 49 : if (mesh->allow_renumbering())
84 27 : mesh->renumber_nodes_and_elements();
85 :
86 : // We need to update the global_boundary_ids, and this is faster
87 : // than a full prepare_for_use()
88 49 : boundary_info.regenerate_id_sets();
89 :
90 98 : return dynamic_pointer_cast<MeshBase>(mesh);
91 49 : }
92 :
93 : BreakMeshByElementGenerator::NodeToElemMapType
94 49 : BreakMeshByElementGenerator::buildSubdomainRestrictedNodeToElemMap(
95 : std::unique_ptr<MeshBase> & mesh, const std::vector<SubdomainID> & subdomains) const
96 : {
97 49 : NodeToElemMapType node_to_elem_map;
98 54167 : for (const auto & elem : mesh->active_element_ptr_range())
99 : {
100 : // Skip if subdomains are specified and the element is not in them
101 28422 : if (!subdomains.empty() &&
102 28422 : std::find(subdomains.begin(), subdomains.end(), elem->subdomain_id()) == subdomains.end())
103 192 : continue;
104 :
105 26867 : std::set<const Elem *> neighbors;
106 26867 : elem->find_point_neighbors(neighbors);
107 :
108 134847 : for (auto n : make_range(elem->n_nodes()))
109 : {
110 : // if ANY neighboring element that contains this node is not in specified subdomains,
111 : // don't add this node to the map, i.e. don't split this node.
112 107980 : bool should_duplicate = true;
113 107980 : if (!subdomains.empty())
114 43168 : for (auto neighbor : neighbors)
115 55533 : if (neighbor->contains_point(elem->node_ref(n)) &&
116 17241 : std::find(subdomains.begin(), subdomains.end(), neighbor->subdomain_id()) ==
117 55533 : subdomains.end())
118 : {
119 320 : should_duplicate = false;
120 320 : break;
121 : }
122 :
123 107980 : if (should_duplicate)
124 107660 : node_to_elem_map[elem->node_id(n)].insert(elem->id());
125 : }
126 26916 : }
127 :
128 49 : return node_to_elem_map;
129 0 : }
130 :
131 : void
132 49 : BreakMeshByElementGenerator::duplicateNodes(std::unique_ptr<MeshBase> & mesh,
133 : const NodeToElemMapType & node_to_elem_map) const
134 : {
135 8740 : for (const auto & [node_id, connected_elem_ids] : node_to_elem_map)
136 116351 : for (auto & connected_elem_id : connected_elem_ids)
137 107660 : if (connected_elem_id != *connected_elem_ids.begin())
138 98969 : duplicateNode(mesh, mesh->elem_ptr(connected_elem_id), mesh->node_ptr(node_id));
139 49 : }
140 :
141 : void
142 98969 : BreakMeshByElementGenerator::duplicateNode(std::unique_ptr<MeshBase> & mesh,
143 : Elem * elem,
144 : const Node * node) const
145 : {
146 98969 : std::unique_ptr<Node> new_node = Node::build(*node, Node::invalid_id);
147 98969 : Node * added_node = mesh->add_node(std::move(new_node));
148 98969 : added_node->processor_id() = elem->processor_id();
149 249021 : for (const auto j : elem->node_index_range())
150 249021 : if (elem->node_id(j) == node->id())
151 : {
152 98969 : elem->set_node(j, added_node);
153 98969 : break;
154 : }
155 :
156 : // Add boundary info to the new node
157 98969 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
158 98969 : std::vector<boundary_id_type> node_boundary_ids;
159 98969 : boundary_info.boundary_ids(node, node_boundary_ids);
160 98969 : boundary_info.add_node(added_node, node_boundary_ids);
161 98969 : }
162 :
163 : void
164 49 : BreakMeshByElementGenerator::createInterface(MeshBase & mesh,
165 : const NodeToElemMapType & node_to_elem_map) const
166 : {
167 49 : std::set<std::pair<dof_id_type, unsigned int>> sides_breaking;
168 :
169 8740 : for (const auto & node_to_elems : node_to_elem_map)
170 116351 : for (const auto & elem_id_i : node_to_elems.second)
171 : {
172 107660 : Elem * elem_i = mesh.elem_ptr(elem_id_i);
173 1986420 : for (const auto & elem_id_j : node_to_elems.second)
174 : {
175 1878760 : Elem * elem_j = mesh.elem_ptr(elem_id_j);
176 1878760 : if (elem_i != elem_j && elem_i->has_neighbor(elem_j))
177 290384 : sides_breaking.insert(std::make_pair(elem_id_i, elem_i->which_neighbor_am_i(elem_j)));
178 : }
179 : }
180 :
181 49 : if (_interface_sides)
182 : {
183 49 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
184 :
185 49 : const auto & existing_boundary_ids = boundary_info.get_boundary_ids();
186 : const boundary_id_type interface_id =
187 49 : existing_boundary_ids.empty() ? 0 : *existing_boundary_ids.rbegin() + 1;
188 49 : boundary_info.sideset_name(interface_id) = _interface_name;
189 :
190 97941 : for (const auto & [elem_id, side] : sides_breaking)
191 97892 : if (_interface_sides > 1 || elem_id > mesh.elem_ptr(elem_id)->neighbor_ptr(side)->id())
192 48946 : boundary_info.add_side(elem_id, side, interface_id);
193 : }
194 :
195 : // Remove element neighbor connections on the broken sides
196 97941 : for (const auto & [elem_id, side] : sides_breaking)
197 97892 : mesh.elem_ref(elem_id).set_neighbor(side, nullptr);
198 49 : }
|
