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 "BreakMeshByBlockGenerator.h"
11 : #include "CastUniquePointer.h"
12 : #include "MooseMeshUtils.h"
13 :
14 : #include "libmesh/distributed_mesh.h"
15 : #include "libmesh/elem.h"
16 : #include "libmesh/partitioner.h"
17 : #include <typeinfo>
18 :
19 : registerMooseObject("MooseApp", BreakMeshByBlockGenerator);
20 :
21 : InputParameters
22 14801 : BreakMeshByBlockGenerator::validParams()
23 : {
24 14801 : InputParameters params = BreakMeshByBlockGeneratorBase::validParams();
25 14801 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
26 14801 : params.addClassDescription(
27 : "Break the mesh at interfaces between blocks. New nodes will be generated so elements on "
28 : "each side of the break are no longer connected. At the moment, this only works on a "
29 : "REPLICATED mesh");
30 14801 : params.addParam<std::vector<SubdomainName>>(
31 : "surrounding_blocks",
32 : "The list of subdomain names surrounding which interfaces will be generated.");
33 14801 : params.addParam<std::vector<std::vector<SubdomainName>>>(
34 : "block_pairs", "The list of subdomain pairs between which interfaces will be generated.");
35 44403 : params.addParam<bool>("add_transition_interface",
36 29602 : false,
37 : "If true and block is not empty, a special boundary named "
38 : "interface_transition is generate between listed blocks and other blocks.");
39 44403 : params.addParam<bool>(
40 29602 : "split_transition_interface", false, "Whether to split the transition interface by blocks.");
41 44403 : params.addParam<bool>("add_interface_on_two_sides",
42 29602 : false,
43 : "Whether to add an additional interface boundary at the other side.");
44 14801 : params.addParam<BoundaryName>(
45 : "interface_transition_name",
46 : "interface_transition",
47 : "the name of the interface transition boundary created when blocks are provided");
48 14801 : return params;
49 0 : }
50 :
51 268 : BreakMeshByBlockGenerator::BreakMeshByBlockGenerator(const InputParameters & parameters)
52 : : BreakMeshByBlockGeneratorBase(parameters),
53 268 : _input(getMesh("input")),
54 268 : _block_pairs_restricted(parameters.isParamSetByUser("block_pairs")),
55 268 : _surrounding_blocks_restricted(parameters.isParamSetByUser("surrounding_blocks")),
56 268 : _add_transition_interface(getParam<bool>("add_transition_interface")),
57 268 : _split_transition_interface(getParam<bool>("split_transition_interface")),
58 268 : _interface_transition_name(getParam<BoundaryName>("interface_transition_name")),
59 804 : _add_interface_on_two_sides(getParam<bool>("add_interface_on_two_sides"))
60 : {
61 268 : if (_block_pairs_restricted && _surrounding_blocks_restricted)
62 0 : paramError("block_pairs_restricted",
63 : "BreakMeshByBlockGenerator: 'surrounding_blocks' and 'block_pairs' can not be used "
64 : "at the same time.");
65 :
66 268 : if (!_add_transition_interface && _split_transition_interface)
67 0 : paramError("split_transition_interface",
68 : "BreakMeshByBlockGenerator cannot split the transition interface because "
69 : "add_transition_interface is false");
70 :
71 268 : if (_add_transition_interface && _block_pairs_restricted)
72 0 : paramError("add_transition_interface",
73 : "BreakMeshByBlockGenerator cannot split the transition interface when 'block_pairs' "
74 : "are specified.");
75 268 : }
76 :
77 : std::unique_ptr<MeshBase>
78 268 : BreakMeshByBlockGenerator::generate()
79 : {
80 268 : std::unique_ptr<MeshBase> mesh = std::move(_input);
81 268 : if (!mesh->is_replicated())
82 0 : mooseError("BreakMeshByBlockGenerator is not implemented for distributed meshes");
83 :
84 : // Try to trick the rest of the world into thinking we're prepared
85 268 : mesh->prepare_for_use();
86 :
87 268 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
88 :
89 : // Handle block restrictions
90 268 : if (_block_pairs_restricted)
91 : {
92 67 : for (const auto & block_name_pair :
93 224 : getParam<std::vector<std::vector<SubdomainName>>>("block_pairs"))
94 : {
95 94 : if (block_name_pair.size() != 2)
96 0 : paramError("block_pairs",
97 : "Each row of 'block_pairs' must have a size of two (block names).");
98 :
99 : // check that the blocks exist in the mesh
100 278 : for (const auto & name : block_name_pair)
101 188 : if (!MooseMeshUtils::hasSubdomainName(*mesh, name))
102 4 : paramError("block_pairs", "The block '", name, "' was not found in the mesh");
103 :
104 90 : const auto block_pair = MooseMeshUtils::getSubdomainIDs(*mesh, block_name_pair);
105 180 : std::pair<SubdomainID, SubdomainID> pair = std::make_pair(
106 180 : std::min(block_pair[0], block_pair[1]), std::max(block_pair[0], block_pair[1]));
107 :
108 90 : _block_pairs.insert(pair);
109 90 : std::copy(block_pair.begin(), block_pair.end(), std::inserter(_block_set, _block_set.end()));
110 90 : }
111 : }
112 201 : else if (_surrounding_blocks_restricted)
113 : {
114 : // check that the blocks exist in the mesh
115 1444 : for (const auto & name : getParam<std::vector<SubdomainName>>("surrounding_blocks"))
116 1399 : if (!MooseMeshUtils::hasSubdomainName(*mesh, name))
117 4 : paramError("surrounding_blocks", "The block '", name, "' was not found in the mesh");
118 :
119 : const auto surrounding_block_ids = MooseMeshUtils::getSubdomainIDs(
120 45 : *mesh, getParam<std::vector<SubdomainName>>("surrounding_blocks"));
121 90 : std::copy(surrounding_block_ids.begin(),
122 : surrounding_block_ids.end(),
123 45 : std::inserter(_block_set, _block_set.end()));
124 45 : }
125 :
126 : // check that a boundary named _interface_transition_name does not already exist in the mesh
127 287 : if ((_block_pairs_restricted || _surrounding_blocks_restricted) && _add_transition_interface &&
128 27 : boundary_info.get_id_by_name(_interface_transition_name) != Moose::INVALID_BOUNDARY_ID)
129 0 : paramError("interface_transition_name",
130 : "BreakMeshByBlockGenerator the specified interface transition boundary name "
131 : "already exits");
132 :
133 : // initialize the node to element map
134 260 : std::map<dof_id_type, std::vector<dof_id_type>> node_to_elem_map;
135 41874 : for (const auto & elem : mesh->active_element_ptr_range())
136 120642 : for (unsigned int n = 0; n < elem->n_nodes(); n++)
137 100095 : node_to_elem_map[elem->node_id(n)].push_back(elem->id());
138 :
139 18429 : for (auto node_it = node_to_elem_map.begin(); node_it != node_to_elem_map.end(); ++node_it)
140 : {
141 18169 : const dof_id_type current_node_id = node_it->first;
142 18169 : const Node * current_node = mesh->node_ptr(current_node_id);
143 :
144 18169 : if (current_node != nullptr)
145 : {
146 : // find node multiplicity
147 18169 : std::set<subdomain_id_type> connected_blocks;
148 118004 : for (auto elem_id = node_it->second.begin(); elem_id != node_it->second.end(); elem_id++)
149 : {
150 99835 : const Elem * current_elem = mesh->elem_ptr(*elem_id);
151 99835 : subdomain_id_type block_id = blockRestrictedElementSubdomainID(current_elem);
152 99835 : if (!_block_pairs_restricted)
153 93067 : connected_blocks.insert(block_id);
154 : else
155 : {
156 6768 : if (block_id != Elem::invalid_subdomain_id)
157 6480 : connected_blocks.insert(block_id);
158 : }
159 : }
160 :
161 18169 : unsigned int node_multiplicity = connected_blocks.size();
162 :
163 : // check if current_node need to be duplicated
164 18169 : if (node_multiplicity > 1)
165 : {
166 : // retrieve connected elements from the map
167 8792 : const std::vector<dof_id_type> & connected_elems = node_it->second;
168 :
169 : // find reference_subdomain_id (e.g. the subdomain with lower id or the
170 : // Elem::invalid_subdomain_id)
171 8792 : auto subdomain_it = connected_blocks.begin();
172 : subdomain_id_type reference_subdomain_id =
173 8792 : connected_blocks.find(Elem::invalid_subdomain_id) != connected_blocks.end()
174 15154 : ? Elem::invalid_subdomain_id
175 6362 : : *subdomain_it;
176 :
177 : // For the block_pairs option, the node that is only shared by specific block pairs will be
178 : // newly created.
179 8792 : bool should_create_new_node = true;
180 8792 : if (_block_pairs_restricted)
181 : {
182 513 : auto elems = node_to_elem_map[current_node->id()];
183 513 : should_create_new_node = false;
184 513 : std::set<subdomain_id_type> sets_blocks_for_this_node;
185 2493 : for (auto elem_id = elems.begin(); elem_id != elems.end(); elem_id++)
186 1980 : sets_blocks_for_this_node.insert(
187 1980 : blockRestrictedElementSubdomainID(mesh->elem_ptr(*elem_id)));
188 513 : if (sets_blocks_for_this_node.size() == 2)
189 : {
190 459 : auto setIt = sets_blocks_for_this_node.begin();
191 459 : if (findBlockPairs(*setIt, *std::next(setIt, 1)))
192 396 : should_create_new_node = true;
193 : }
194 513 : }
195 :
196 : // multiplicity counter to keep track of how many nodes we added
197 8792 : unsigned int multiplicity_counter = node_multiplicity;
198 68289 : for (auto elem_id : connected_elems)
199 : {
200 : // all the duplicate nodes are added and assigned
201 59506 : if (multiplicity_counter == 0)
202 9 : break;
203 :
204 59497 : Elem * current_elem = mesh->elem_ptr(elem_id);
205 59497 : subdomain_id_type block_id = blockRestrictedElementSubdomainID(current_elem);
206 :
207 85299 : if ((blockRestrictedElementSubdomainID(current_elem) != reference_subdomain_id) ||
208 25802 : (_block_pairs_restricted & findBlockPairs(reference_subdomain_id, block_id)))
209 : {
210 : // assign the newly added node to current_elem
211 33695 : Node * new_node = nullptr;
212 :
213 33695 : std::vector<boundary_id_type> node_boundary_ids;
214 :
215 160937 : for (unsigned int node_id = 0; node_id < current_elem->n_nodes(); ++node_id)
216 144701 : if (current_elem->node_id(node_id) ==
217 144701 : current_node->id()) // if current node == node on element
218 : {
219 17459 : if (should_create_new_node)
220 : {
221 : // add new node
222 17279 : new_node = Node::build(*current_node, mesh->n_nodes()).release();
223 :
224 : // We're duplicating nodes so that each subdomain elem has its own copy, so it
225 : // seems natural to assign this new node the same proc id as corresponding
226 : // subdomain elem
227 17279 : new_node->processor_id() = current_elem->processor_id();
228 17279 : mesh->add_node(new_node);
229 17279 : current_elem->set_node(node_id, new_node);
230 : // Add boundary info to the new node
231 17279 : boundary_info.boundary_ids(current_node, node_boundary_ids);
232 17279 : boundary_info.add_node(new_node, node_boundary_ids);
233 : }
234 :
235 17459 : multiplicity_counter--; // node created, update multiplicity counter
236 :
237 17459 : break; // ones the proper node has been fixed in one element we can break the
238 : // loop
239 : }
240 :
241 33695 : if (should_create_new_node)
242 : {
243 473685 : for (auto connected_elem_id : connected_elems)
244 : {
245 440170 : Elem * connected_elem = mesh->elem_ptr(connected_elem_id);
246 :
247 : // Assign the newly added node to other connected elements with the same
248 : // block_id
249 440170 : if (connected_elem->subdomain_id() == current_elem->subdomain_id() &&
250 : connected_elem != current_elem)
251 : {
252 714780 : for (unsigned int node_id = 0; node_id < connected_elem->n_nodes(); ++node_id)
253 582318 : if (connected_elem->node_id(node_id) ==
254 582318 : current_node->id()) // if current node == node on element
255 : {
256 16236 : connected_elem->set_node(node_id, new_node);
257 16236 : break;
258 : }
259 : }
260 : }
261 : }
262 33695 : }
263 : }
264 :
265 : // create blocks pair and assign element side to new interface boundary map
266 68307 : for (auto elem_id : connected_elems)
267 : {
268 837912 : for (auto connected_elem_id : connected_elems)
269 : {
270 778397 : Elem * current_elem = mesh->elem_ptr(elem_id);
271 778397 : Elem * connected_elem = mesh->elem_ptr(connected_elem_id);
272 778397 : subdomain_id_type curr_elem_subid = blockRestrictedElementSubdomainID(current_elem);
273 : subdomain_id_type connected_elem_subid =
274 778397 : blockRestrictedElementSubdomainID(connected_elem);
275 :
276 778397 : if (current_elem != connected_elem && curr_elem_subid < connected_elem_subid)
277 : {
278 218609 : if (current_elem->has_neighbor(connected_elem))
279 : {
280 34181 : dof_id_type elem_id = current_elem->id();
281 34181 : dof_id_type connected_elem_id = connected_elem->id();
282 34181 : unsigned int side = current_elem->which_neighbor_am_i(connected_elem);
283 : unsigned int connected_elem_side =
284 34181 : connected_elem->which_neighbor_am_i(current_elem);
285 :
286 : // in this case we need to play a game to reorder the sides
287 34181 : if (_block_pairs_restricted || _surrounding_blocks_restricted)
288 : {
289 18144 : connected_elem_subid = connected_elem->subdomain_id();
290 18144 : if (curr_elem_subid > connected_elem_subid) // we need to switch the ids
291 : {
292 3258 : connected_elem_subid = current_elem->subdomain_id();
293 3258 : curr_elem_subid = connected_elem->subdomain_id();
294 :
295 3258 : elem_id = connected_elem->id();
296 3258 : side = connected_elem->which_neighbor_am_i(current_elem);
297 :
298 3258 : connected_elem_id = current_elem->id();
299 3258 : connected_elem_side = current_elem->which_neighbor_am_i(connected_elem);
300 : }
301 : }
302 :
303 : std::pair<subdomain_id_type, subdomain_id_type> blocks_pair =
304 34181 : std::make_pair(curr_elem_subid, connected_elem_subid);
305 :
306 : std::pair<subdomain_id_type, subdomain_id_type> blocks_pair2 =
307 34181 : std::make_pair(connected_elem_subid, curr_elem_subid);
308 :
309 34181 : if (_block_pairs_restricted)
310 : {
311 1044 : if (findBlockPairs(blockRestrictedElementSubdomainID(current_elem),
312 1044 : blockRestrictedElementSubdomainID(connected_elem)))
313 : {
314 864 : _new_boundary_sides_map[blocks_pair].insert(std::make_pair(elem_id, side));
315 864 : if (_add_interface_on_two_sides)
316 1368 : _new_boundary_sides_map[blocks_pair2].insert(
317 1368 : std::make_pair(connected_elem_id, connected_elem_side));
318 : }
319 : }
320 : else
321 : {
322 33137 : _new_boundary_sides_map[blocks_pair].insert(std::make_pair(elem_id, side));
323 33137 : if (_add_interface_on_two_sides)
324 0 : _new_boundary_sides_map[blocks_pair2].insert(
325 0 : std::make_pair(connected_elem_id, connected_elem_side));
326 : }
327 : }
328 : }
329 : }
330 : }
331 :
332 : } // end multiplicity check
333 18169 : } // end loop over nodes
334 : } // end nodeptr check
335 :
336 260 : addInterfaceBoundary(*mesh);
337 260 : Partitioner::set_node_processor_ids(*mesh);
338 520 : return dynamic_pointer_cast<MeshBase>(mesh);
339 260 : }
340 :
341 : void
342 260 : BreakMeshByBlockGenerator::addInterfaceBoundary(MeshBase & mesh)
343 : {
344 260 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
345 :
346 : boundary_id_type boundary_id;
347 260 : boundary_id_type boundary_id_interface = Moose::INVALID_BOUNDARY_ID;
348 260 : boundary_id_type boundary_id_interface_transition = Moose::INVALID_BOUNDARY_ID;
349 :
350 260 : BoundaryName boundary_name;
351 :
352 : // loop over boundary sides
353 5659 : for (auto & boundary_side_map : _new_boundary_sides_map)
354 : {
355 5399 : if (!(_block_pairs_restricted || _surrounding_blocks_restricted) ||
356 4401 : ((_block_pairs_restricted || _surrounding_blocks_restricted) && !_add_transition_interface))
357 : {
358 : // find the appropriate boundary name and id
359 : // given primary and secondary block ID
360 2834 : if (_split_interface)
361 1485 : findBoundaryNameAndInd(mesh,
362 1485 : boundary_side_map.first.first,
363 1485 : boundary_side_map.first.second,
364 : boundary_name,
365 : boundary_id,
366 : boundary_info);
367 : else
368 : {
369 1349 : boundary_name = _interface_name;
370 2698 : boundary_id_interface = boundary_id_interface == Moose::INVALID_BOUNDARY_ID
371 1349 : ? findFreeBoundaryId(mesh)
372 : : boundary_id_interface;
373 1349 : boundary_id = boundary_id_interface;
374 1349 : boundary_info.sideset_name(boundary_id_interface) = boundary_name;
375 : }
376 : }
377 : else // block resticted with transition boundary
378 : {
379 :
380 : const bool interior_boundary =
381 4644 : _block_set.find(boundary_side_map.first.first) != _block_set.end() &&
382 2079 : _block_set.find(boundary_side_map.first.second) != _block_set.end();
383 :
384 2565 : if ((_split_interface && interior_boundary) ||
385 1485 : (_split_transition_interface && !interior_boundary))
386 : {
387 1710 : findBoundaryNameAndInd(mesh,
388 1710 : boundary_side_map.first.first,
389 1710 : boundary_side_map.first.second,
390 : boundary_name,
391 : boundary_id,
392 : boundary_info);
393 : }
394 855 : else if (interior_boundary)
395 : {
396 540 : boundary_name = _interface_name;
397 1080 : boundary_id_interface = boundary_id_interface == Moose::INVALID_BOUNDARY_ID
398 540 : ? findFreeBoundaryId(mesh)
399 : : boundary_id_interface;
400 :
401 540 : boundary_id = boundary_id_interface;
402 540 : boundary_info.sideset_name(boundary_id_interface) = boundary_name;
403 : }
404 : else
405 : {
406 315 : boundary_id_interface_transition =
407 315 : boundary_id_interface_transition == Moose::INVALID_BOUNDARY_ID
408 315 : ? findFreeBoundaryId(mesh)
409 : : boundary_id_interface_transition;
410 315 : boundary_id = boundary_id_interface_transition;
411 315 : boundary_info.sideset_name(boundary_id) = _interface_transition_name;
412 : }
413 : }
414 : // loop over all the side belonging to each pair and add it to the proper interface
415 15190 : for (auto & element_side : boundary_side_map.second)
416 9791 : boundary_info.add_side(element_side.first, element_side.second, boundary_id);
417 : }
418 260 : }
419 :
420 : subdomain_id_type
421 1779691 : BreakMeshByBlockGenerator::blockRestrictedElementSubdomainID(const Elem * elem)
422 : {
423 1779691 : subdomain_id_type elem_subdomain_id = elem->subdomain_id();
424 2083171 : if ((_block_pairs_restricted || _surrounding_blocks_restricted) &&
425 2083171 : (_block_set.find(elem_subdomain_id) == _block_set.end()))
426 99756 : elem_subdomain_id = Elem::invalid_subdomain_id;
427 :
428 1779691 : return elem_subdomain_id;
429 : }
430 :
431 : bool
432 27305 : BreakMeshByBlockGenerator::findBlockPairs(subdomain_id_type block_one, subdomain_id_type block_two)
433 : {
434 29060 : for (auto block_pair : _block_pairs)
435 3015 : if ((block_pair.first == block_one && block_pair.second == block_two) ||
436 1755 : (block_pair.first == block_two && block_pair.second == block_one))
437 1260 : return true;
438 26045 : return false;
439 : }
|
