www.mooseframework.org
BreakMeshByBlockGenerator.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 
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 
20 
23 {
25  params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
26  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  params.addParam<std::vector<SubdomainName>>(
31  "surrounding_blocks",
32  "The list of subdomain names surrounding which interfaces will be generated.");
33  params.addParam<std::vector<std::vector<SubdomainName>>>(
34  "block_pairs", "The list of subdomain pairs between which interfaces will be generated.");
35  params.addParam<bool>("add_transition_interface",
36  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  params.addParam<bool>(
40  "split_transition_interface", false, "Whether to split the transition interface by blocks.");
41  params.addParam<bool>("add_interface_on_two_sides",
42  false,
43  "Whether to add an additional interface boundary at the other side.");
44  params.addParam<BoundaryName>(
45  "interface_transition_name",
46  "interface_transition",
47  "the name of the interface transition boundary created when blocks are provided");
48  return params;
49 }
50 
52  : BreakMeshByBlockGeneratorBase(parameters),
53  _input(getMesh("input")),
54  _block_pairs_restricted(parameters.isParamSetByUser("block_pairs")),
55  _surrounding_blocks_restricted(parameters.isParamSetByUser("surrounding_blocks")),
56  _add_transition_interface(getParam<bool>("add_transition_interface")),
57  _split_transition_interface(getParam<bool>("split_transition_interface")),
58  _interface_transition_name(getParam<BoundaryName>("interface_transition_name")),
59  _add_interface_on_two_sides(getParam<bool>("add_interface_on_two_sides"))
60 {
62  paramError("block_pairs_restricted",
63  "BreakMeshByBlockGenerator: 'surrounding_blocks' and 'block_pairs' can not be used "
64  "at the same time.");
65 
67  paramError("split_transition_interface",
68  "BreakMeshByBlockGenerator cannot split the transition interface because "
69  "add_transition_interface is false");
70 
72  paramError("add_transition_interface",
73  "BreakMeshByBlockGenerator cannot split the transition interface when 'block_pairs' "
74  "are specified.");
75 }
76 
77 std::unique_ptr<MeshBase>
79 {
80  std::unique_ptr<MeshBase> mesh = std::move(_input);
81  if (!mesh->is_replicated())
82  mooseError("BreakMeshByBlockGenerator is not implemented for distributed meshes");
83 
84  // Try to trick the rest of the world into thinking we're prepared
85  mesh->prepare_for_use();
86 
87  BoundaryInfo & boundary_info = mesh->get_boundary_info();
88 
89  // Handle block restrictions
91  {
92  for (const auto & block_name_pair :
93  getParam<std::vector<std::vector<SubdomainName>>>("block_pairs"))
94  {
95  if (block_name_pair.size() != 2)
96  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  for (const auto & name : block_name_pair)
102  paramError("block_pairs", "The block '", name, "' was not found in the mesh");
103 
104  const auto block_pair = MooseMeshUtils::getSubdomainIDs(*mesh, block_name_pair);
105  std::pair<SubdomainID, SubdomainID> pair = std::make_pair(
106  std::min(block_pair[0], block_pair[1]), std::max(block_pair[0], block_pair[1]));
107 
108  _block_pairs.insert(pair);
109  std::copy(block_pair.begin(), block_pair.end(), std::inserter(_block_set, _block_set.end()));
110  }
111  }
113  {
114  // check that the blocks exist in the mesh
115  for (const auto & name : getParam<std::vector<SubdomainName>>("surrounding_blocks"))
117  paramError("surrounding_blocks", "The block '", name, "' was not found in the mesh");
118 
119  const auto surrounding_block_ids = MooseMeshUtils::getSubdomainIDs(
120  *mesh, getParam<std::vector<SubdomainName>>("surrounding_blocks"));
121  std::copy(surrounding_block_ids.begin(),
122  surrounding_block_ids.end(),
123  std::inserter(_block_set, _block_set.end()));
124  }
125 
126  // check that a boundary named _interface_transition_name does not already exist in the mesh
128  boundary_info.get_id_by_name(_interface_transition_name) != Moose::INVALID_BOUNDARY_ID)
129  paramError("interface_transition_name",
130  "BreakMeshByBlockGenerator the specified interface transition boundary name "
131  "already exits");
132 
133  // initialize the node to element map
134  std::map<dof_id_type, std::vector<dof_id_type>> node_to_elem_map;
135  for (const auto & elem : mesh->active_element_ptr_range())
136  for (unsigned int n = 0; n < elem->n_nodes(); n++)
137  node_to_elem_map[elem->node_id(n)].push_back(elem->id());
138 
139  for (auto node_it = node_to_elem_map.begin(); node_it != node_to_elem_map.end(); ++node_it)
140  {
141  const dof_id_type current_node_id = node_it->first;
142  const Node * current_node = mesh->node_ptr(current_node_id);
143 
144  if (current_node != nullptr)
145  {
146  // find node multiplicity
147  std::set<subdomain_id_type> connected_blocks;
148  for (auto elem_id = node_it->second.begin(); elem_id != node_it->second.end(); elem_id++)
149  {
150  const Elem * current_elem = mesh->elem_ptr(*elem_id);
151  subdomain_id_type block_id = blockRestrictedElementSubdomainID(current_elem);
153  connected_blocks.insert(block_id);
154  else
155  {
156  if (block_id != Elem::invalid_subdomain_id)
157  connected_blocks.insert(block_id);
158  }
159  }
160 
161  unsigned int node_multiplicity = connected_blocks.size();
162 
163  // check if current_node need to be duplicated
164  if (node_multiplicity > 1)
165  {
166  // retrieve connected elements from the map
167  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  auto subdomain_it = connected_blocks.begin();
172  subdomain_id_type reference_subdomain_id =
173  connected_blocks.find(Elem::invalid_subdomain_id) != connected_blocks.end()
174  ? Elem::invalid_subdomain_id
175  : *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  bool should_create_new_node = true;
181  {
182  auto elems = node_to_elem_map[current_node->id()];
183  should_create_new_node = false;
184  std::set<subdomain_id_type> sets_blocks_for_this_node;
185  for (auto elem_id = elems.begin(); elem_id != elems.end(); elem_id++)
186  sets_blocks_for_this_node.insert(
187  blockRestrictedElementSubdomainID(mesh->elem_ptr(*elem_id)));
188  if (sets_blocks_for_this_node.size() == 2)
189  {
190  auto setIt = sets_blocks_for_this_node.begin();
191  if (findBlockPairs(*setIt, *std::next(setIt, 1)))
192  should_create_new_node = true;
193  }
194  }
195 
196  // multiplicity counter to keep track of how many nodes we added
197  unsigned int multiplicity_counter = node_multiplicity;
198  for (auto elem_id : connected_elems)
199  {
200  // all the duplicate nodes are added and assigned
201  if (multiplicity_counter == 0)
202  break;
203 
204  Elem * current_elem = mesh->elem_ptr(elem_id);
205  subdomain_id_type block_id = blockRestrictedElementSubdomainID(current_elem);
206 
207  if ((blockRestrictedElementSubdomainID(current_elem) != reference_subdomain_id) ||
208  (_block_pairs_restricted & findBlockPairs(reference_subdomain_id, block_id)))
209  {
210  // assign the newly added node to current_elem
211  Node * new_node = nullptr;
212 
213  std::vector<boundary_id_type> node_boundary_ids;
214 
215  for (unsigned int node_id = 0; node_id < current_elem->n_nodes(); ++node_id)
216  if (current_elem->node_id(node_id) ==
217  current_node->id()) // if current node == node on element
218  {
219  if (should_create_new_node)
220  {
221  // add new node
222  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  new_node->processor_id() = current_elem->processor_id();
228  mesh->add_node(new_node);
229  current_elem->set_node(node_id) = new_node;
230  // Add boundary info to the new node
231  boundary_info.boundary_ids(current_node, node_boundary_ids);
232  boundary_info.add_node(new_node, node_boundary_ids);
233  }
234 
235  multiplicity_counter--; // node created, update multiplicity counter
236 
237  break; // ones the proper node has been fixed in one element we can break the
238  // loop
239  }
240 
241  if (should_create_new_node)
242  {
243  for (auto connected_elem_id : connected_elems)
244  {
245  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  if (connected_elem->subdomain_id() == current_elem->subdomain_id() &&
250  connected_elem != current_elem)
251  {
252  for (unsigned int node_id = 0; node_id < connected_elem->n_nodes(); ++node_id)
253  if (connected_elem->node_id(node_id) ==
254  current_node->id()) // if current node == node on element
255  {
256  connected_elem->set_node(node_id) = new_node;
257  break;
258  }
259  }
260  }
261  }
262  }
263  }
264 
265  // create blocks pair and assign element side to new interface boundary map
266  for (auto elem_id : connected_elems)
267  {
268  for (auto connected_elem_id : connected_elems)
269  {
270  Elem * current_elem = mesh->elem_ptr(elem_id);
271  Elem * connected_elem = mesh->elem_ptr(connected_elem_id);
272  subdomain_id_type curr_elem_subid = blockRestrictedElementSubdomainID(current_elem);
273  subdomain_id_type connected_elem_subid =
274  blockRestrictedElementSubdomainID(connected_elem);
275 
276  if (current_elem != connected_elem && curr_elem_subid < connected_elem_subid)
277  {
278  if (current_elem->has_neighbor(connected_elem))
279  {
280  dof_id_type elem_id = current_elem->id();
281  dof_id_type connected_elem_id = connected_elem->id();
282  unsigned int side = current_elem->which_neighbor_am_i(connected_elem);
283  unsigned int connected_elem_side =
284  connected_elem->which_neighbor_am_i(current_elem);
285 
286  // in this case we need to play a game to reorder the sides
288  {
289  connected_elem_subid = connected_elem->subdomain_id();
290  if (curr_elem_subid > connected_elem_subid) // we need to switch the ids
291  {
292  connected_elem_subid = current_elem->subdomain_id();
293  curr_elem_subid = connected_elem->subdomain_id();
294 
295  elem_id = connected_elem->id();
296  side = connected_elem->which_neighbor_am_i(current_elem);
297 
298  connected_elem_id = current_elem->id();
299  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  std::make_pair(curr_elem_subid, connected_elem_subid);
305 
306  std::pair<subdomain_id_type, subdomain_id_type> blocks_pair2 =
307  std::make_pair(connected_elem_subid, curr_elem_subid);
308 
310  {
312  blockRestrictedElementSubdomainID(connected_elem)))
313  {
314  _new_boundary_sides_map[blocks_pair].insert(std::make_pair(elem_id, side));
316  _new_boundary_sides_map[blocks_pair2].insert(
317  std::make_pair(connected_elem_id, connected_elem_side));
318  }
319  }
320  else
321  {
322  _new_boundary_sides_map[blocks_pair].insert(std::make_pair(elem_id, side));
324  _new_boundary_sides_map[blocks_pair2].insert(
325  std::make_pair(connected_elem_id, connected_elem_side));
326  }
327  }
328  }
329  }
330  }
331 
332  } // end multiplicity check
333  } // end loop over nodes
334  } // end nodeptr check
335 
337  Partitioner::set_node_processor_ids(*mesh);
338  return dynamic_pointer_cast<MeshBase>(mesh);
339 }
340 
341 void
343 {
344  BoundaryInfo & boundary_info = mesh.get_boundary_info();
345 
346  boundary_id_type boundary_id;
347  boundary_id_type boundary_id_interface = Moose::INVALID_BOUNDARY_ID;
348  boundary_id_type boundary_id_interface_transition = Moose::INVALID_BOUNDARY_ID;
349 
350  BoundaryName boundary_name;
351 
352  // loop over boundary sides
353  for (auto & boundary_side_map : _new_boundary_sides_map)
354  {
357  {
358  // find the appropriate boundary name and id
359  // given primary and secondary block ID
360  if (_split_interface)
362  boundary_side_map.first.first,
363  boundary_side_map.first.second,
364  boundary_name,
365  boundary_id,
366  boundary_info);
367  else
368  {
369  boundary_name = _interface_name;
370  boundary_id_interface = boundary_id_interface == Moose::INVALID_BOUNDARY_ID
372  : boundary_id_interface;
373  boundary_id = boundary_id_interface;
374  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  _block_set.find(boundary_side_map.first.first) != _block_set.end() &&
382  _block_set.find(boundary_side_map.first.second) != _block_set.end();
383 
384  if ((_split_interface && interior_boundary) ||
385  (_split_transition_interface && !interior_boundary))
386  {
388  boundary_side_map.first.first,
389  boundary_side_map.first.second,
390  boundary_name,
391  boundary_id,
392  boundary_info);
393  }
394  else if (interior_boundary)
395  {
396  boundary_name = _interface_name;
397  boundary_id_interface = boundary_id_interface == Moose::INVALID_BOUNDARY_ID
399  : boundary_id_interface;
400 
401  boundary_id = boundary_id_interface;
402  boundary_info.sideset_name(boundary_id_interface) = boundary_name;
403  }
404  else
405  {
406  boundary_id_interface_transition =
407  boundary_id_interface_transition == Moose::INVALID_BOUNDARY_ID
409  : boundary_id_interface_transition;
410  boundary_id = boundary_id_interface_transition;
411  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  for (auto & element_side : boundary_side_map.second)
416  boundary_info.add_side(element_side.first, element_side.second, boundary_id);
417  }
418 }
419 
422 {
423  subdomain_id_type elem_subdomain_id = elem->subdomain_id();
425  (_block_set.find(elem_subdomain_id) == _block_set.end()))
426  elem_subdomain_id = Elem::invalid_subdomain_id;
427 
428  return elem_subdomain_id;
429 }
430 
431 bool
433 {
434  for (auto block_pair : _block_pairs)
435  if ((block_pair.first == block_one && block_pair.second == block_two) ||
436  (block_pair.first == block_two && block_pair.second == block_one))
437  return true;
438  return false;
439 }
registerMooseObject("MooseApp", BreakMeshByBlockGenerator)
void addInterfaceBoundary(MeshBase &mesh)
generate the new boundary interface
std::unordered_set< std::pair< SubdomainID, SubdomainID > > _block_pairs
set of subdomain pairs between which interfaces will be generated.
static InputParameters validParams()
void findBoundaryNameAndInd(MeshBase &mesh, const subdomain_id_type &, const subdomain_id_type &, std::string &, boundary_id_type &, BoundaryInfo &)
given the primary and secondary blocks this method return the appropriate boundary id and name ...
const BoundaryID INVALID_BOUNDARY_ID
Definition: MooseTypes.C:24
const bool _surrounding_blocks_restricted
whether interfaces will be generated surrounding blocks
const bool _split_transition_interface
whether to split the transition boundary between the blocks and the rest of the mesh ...
TestClass subdomain_id_type
MeshBase & mesh
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
std::unique_ptr< T_DEST, T_DELETER > dynamic_pointer_cast(std::unique_ptr< T_SRC, T_DELETER > &src)
These are reworked from https://stackoverflow.com/a/11003103.
std::vector< subdomain_id_type > getSubdomainIDs(const libMesh::MeshBase &mesh, const std::vector< SubdomainName > &subdomain_name)
Get the associated subdomainIDs for the subdomain names that are passed in.
virtual const std::string & name() const
Get the name of the class.
Definition: MooseBase.h:56
const bool _block_pairs_restricted
whether interfaces will be generated between block pairs
void addRequiredParam(const std::string &name, const std::string &doc_string)
This method adds a parameter and documentation string to the InputParameters object that will be extr...
bool findBlockPairs(subdomain_id_type block_one, subdomain_id_type block_two)
Return true if block_one and block_two are found in users&#39; provided block_pairs list.
auto max(const L &left, const R &right)
int8_t boundary_id_type
const BoundaryName _interface_transition_name
the name of the transition interface
std::string _interface_name
the name of the new interface
std::unique_ptr< MeshBase > & _input
the mesh to modify
const T & getParam(const std::string &name) const
Retrieve a parameter for the object.
bool _split_interface
the flag to split the interface by block
bool hasSubdomainName(MeshBase &input_mesh, const SubdomainName &name)
Whether a particular subdomain name exists in the mesh.
void paramError(const std::string &param, Args... args) const
Emits an error prefixed with the file and line number of the given param (from the input file) along ...
subdomain_id_type blockRestrictedElementSubdomainID(const Elem *elem)
This is a helper method to avoid recoding the same if everywhere.
std::unordered_set< SubdomainID > _block_set
set of the blocks to split the mesh on
const bool _add_transition_interface
whether to add a boundary when splitting the mesh
BreakMeshByBlockGenerator(const InputParameters &parameters)
void mooseError(Args &&... args) const
Emits an error prefixed with object name and type.
const bool _add_interface_on_two_sides
whether to add two sides interface boundaries
void addClassDescription(const std::string &doc_string)
This method adds a description of the class that will be displayed in the input file syntax dump...
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an option parameter and a documentation string to the InputParameters object...
std::map< std::pair< subdomain_id_type, subdomain_id_type >, std::set< std::pair< dof_id_type, unsigned int > > > _new_boundary_sides_map
auto min(const L &left, const R &right)
BoundaryID findFreeBoundaryId(MeshBase &mesh)
this method finds the first free boundary id
std::unique_ptr< MeshBase > generate() override
Generate / modify the mesh.
uint8_t dof_id_type