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 "CombinerGenerator.h"
11 :
12 : #include "CastUniquePointer.h"
13 : #include "MooseUtils.h"
14 : #include "DelimitedFileReader.h"
15 :
16 : #include "libmesh/replicated_mesh.h"
17 : #include "libmesh/unstructured_mesh.h"
18 : #include "libmesh/mesh_modification.h"
19 : #include "libmesh/point.h"
20 : #include "libmesh/node.h"
21 :
22 : registerMooseObject("MooseApp", CombinerGenerator);
23 :
24 : InputParameters
25 15669 : CombinerGenerator::validParams()
26 : {
27 15669 : InputParameters params = MeshGenerator::validParams();
28 :
29 15669 : params.addClassDescription(
30 : "Combine multiple meshes (or copies of one mesh) together into one (disjoint) mesh. Can "
31 : "optionally translate those meshes before combining them.");
32 :
33 15669 : params.addRequiredParam<std::vector<MeshGeneratorName>>(
34 : "inputs",
35 : "The input MeshGenerators. This can either be N generators or 1 generator. If only 1 is "
36 : "given then 'positions' must also be given.");
37 :
38 15669 : params.addParam<std::vector<Point>>(
39 : "positions",
40 : "The (optional) position of each given mesh. If N 'inputs' were given then this must either "
41 : "be left blank or N positions must be given. If 1 input was given then this MUST be "
42 : "provided.");
43 :
44 15669 : params.addParam<std::vector<FileName>>(
45 : "positions_file", "Alternative way to provide the position of each given mesh.");
46 :
47 47007 : params.addParam<bool>("avoid_merging_subdomains",
48 31338 : false,
49 : "Whether to prevent merging subdomains by offsetting ids. The first mesh "
50 : "in the input will keep the same subdomains ids, the others will have "
51 : "offsets. All subdomain names will remain valid");
52 47007 : params.addParam<bool>("avoid_merging_boundaries",
53 31338 : false,
54 : "Whether to prevent merging sidesets by offsetting ids. The first mesh "
55 : "in the input will keep the same boundary ids, the others will have "
56 : "offsets. All boundary names will remain valid");
57 :
58 15669 : return params;
59 0 : }
60 :
61 702 : CombinerGenerator::CombinerGenerator(const InputParameters & parameters)
62 : : MeshGenerator(parameters),
63 702 : _meshes(getMeshes("inputs")),
64 702 : _input_names(getParam<std::vector<MeshGeneratorName>>("inputs")),
65 702 : _avoid_merging_subdomains(getParam<bool>("avoid_merging_subdomains")),
66 2106 : _avoid_merging_boundaries(getParam<bool>("avoid_merging_boundaries"))
67 : {
68 702 : if (_input_names.empty())
69 0 : paramError("input_names", "You need to specify at least one MeshGenerator as an input.");
70 :
71 702 : if (isParamValid("positions") && isParamValid("positions_file"))
72 0 : mooseError("Both 'positions' and 'positions_file' cannot be specified simultaneously in "
73 : "CombinerGenerator ",
74 0 : _name);
75 :
76 702 : if (_input_names.size() == 1)
77 51 : if (!isParamValid("positions") && !isParamValid("positions_file"))
78 0 : paramError("positions",
79 : "If only one input mesh is given, then 'positions' or 'positions_file' must also "
80 : "be supplied");
81 702 : }
82 :
83 : void
84 683 : CombinerGenerator::fillPositions()
85 : {
86 683 : if (isParamValid("positions"))
87 : {
88 264 : _positions = getParam<std::vector<Point>>("positions");
89 :
90 : // the check in the constructor wont catch error where the user sets positions = ''
91 264 : if ((_input_names.size() == 1) && _positions.empty())
92 4 : paramError("positions", "If only one input mesh is given, then 'positions' cannot be empty.");
93 :
94 260 : if (_input_names.size() != 1)
95 229 : if (_positions.size() && (_input_names.size() != _positions.size()))
96 4 : paramError(
97 : "positions",
98 : "If more than one input mesh is provided then the number of positions provided must "
99 : "exactly match the number of input meshes.");
100 : }
101 419 : else if (isParamValid("positions_file"))
102 : {
103 30 : std::vector<FileName> positions_file = getParam<std::vector<FileName>>("positions_file");
104 :
105 : // the check in the constructor wont catch error where the user sets positions_file = ''
106 30 : if ((_input_names.size() == 1) && positions_file.empty())
107 4 : paramError("positions_file",
108 : "If only one input mesh is given, then 'positions_file' cannot be empty.");
109 :
110 48 : for (const auto & f : positions_file)
111 : {
112 26 : MooseUtils::DelimitedFileReader file(f, &_communicator);
113 26 : file.setFormatFlag(MooseUtils::DelimitedFileReader::FormatFlag::ROWS);
114 26 : file.read();
115 :
116 26 : const std::vector<Point> & data = file.getDataAsPoints();
117 :
118 26 : if (_input_names.size() != 1)
119 15 : if (data.size() && (_input_names.size() != data.size()))
120 4 : paramError("positions_file",
121 : "If more than one input mesh is provided then the number of positions must "
122 : "exactly match the number of input meshes.");
123 :
124 88 : for (const auto & d : data)
125 66 : _positions.push_back(d);
126 22 : }
127 22 : }
128 667 : }
129 :
130 : std::unique_ptr<MeshBase>
131 683 : CombinerGenerator::generate()
132 : {
133 : // Two cases:
134 : // 1. Multiple input meshes and optional positions
135 : // 2. One input mesh and multiple positions
136 683 : fillPositions();
137 :
138 : // Case 1
139 667 : if (_meshes.size() != 1)
140 : {
141 : // merge all meshes into the first one
142 625 : auto mesh = dynamic_pointer_cast<UnstructuredMesh>(*_meshes[0]);
143 :
144 625 : if (!mesh)
145 0 : paramError("inputs", _input_names[0], " is not a valid unstructured mesh");
146 :
147 : // Move the first input mesh if applicable
148 625 : if (_positions.size())
149 : {
150 236 : MeshTools::Modification::translate(
151 236 : *mesh, _positions[0](0), _positions[0](1), _positions[0](2));
152 : }
153 :
154 : // Read in all of the other meshes
155 2715 : for (MooseIndex(_meshes) i = 1; i < _meshes.size(); ++i)
156 : {
157 2090 : auto other_mesh = dynamic_pointer_cast<UnstructuredMesh>(*_meshes[i]);
158 :
159 2090 : if (!other_mesh)
160 0 : paramError("inputs", _input_names[i], " is not a valid unstructured mesh");
161 :
162 : // Move It
163 2090 : if (_positions.size())
164 : {
165 774 : MeshTools::Modification::translate(
166 774 : *other_mesh, _positions[i](0), _positions[i](1), _positions[i](2));
167 : }
168 :
169 2090 : copyIntoMesh(*mesh, *other_mesh);
170 2090 : }
171 :
172 625 : mesh->set_isnt_prepared();
173 625 : return dynamic_pointer_cast<MeshBase>(mesh);
174 625 : }
175 : else // Case 2
176 : {
177 42 : auto input_mesh = dynamic_pointer_cast<UnstructuredMesh>(*_meshes[0]);
178 :
179 42 : if (!input_mesh)
180 0 : paramError("inputs", _input_names[0], " is not a valid unstructured mesh");
181 :
182 : // Make a copy and displace it in order to get the final mesh started
183 : auto copy =
184 42 : input_mesh->clone(); // This is required because dynamic_pointer_cast() requires an l-value
185 42 : auto final_mesh = dynamic_pointer_cast<UnstructuredMesh>(copy);
186 :
187 42 : if (!final_mesh)
188 0 : mooseError("Unable to copy mesh!");
189 :
190 42 : MeshTools::Modification::translate(
191 42 : *final_mesh, _positions[0](0), _positions[0](1), _positions[0](2));
192 :
193 : // Here's the way this is going to work:
194 : // I'm going to make one more copy of the input_mesh so that I can move it and copy it in
195 : // Then, after it's copied in I'm going to reset its coordinates by looping over the input_mesh
196 : // and resetting the nodal positions.
197 : // This could be done without the copy - you would translate the mesh then translate it back...
198 : // However, I'm worried about floating point roundoff. If you were doing this 100,000 times or
199 : // more then the mesh could "drift" away from its original position. I really want the
200 : // translations to be exact each time.
201 : // I suppose that it is technically possible to just save off a datastructure (map, etc.) that
202 : // could hold the nodal positions only (instead of a copy of the mesh) but I'm not sure that
203 : // would really save much... we'll see if it shows up in profiling somewhere
204 42 : copy = input_mesh->clone();
205 42 : auto translated_mesh = dynamic_pointer_cast<UnstructuredMesh>(copy);
206 :
207 42 : if (!translated_mesh)
208 0 : mooseError("Unable to copy mesh!");
209 :
210 86 : for (MooseIndex(_meshes) i = 1; i < _positions.size(); ++i)
211 : {
212 : // Move
213 44 : MeshTools::Modification::translate(
214 44 : *translated_mesh, _positions[i](0), _positions[i](1), _positions[i](2));
215 :
216 : // Copy into final mesh
217 44 : copyIntoMesh(*final_mesh, *translated_mesh);
218 :
219 : // Reset nodal coordinates
220 9988 : for (auto translated_node_ptr : translated_mesh->node_ptr_range())
221 : {
222 4972 : auto & translated_node = *translated_node_ptr;
223 4972 : auto & input_node = input_mesh->node_ref(translated_node_ptr->id());
224 :
225 19888 : for (MooseIndex(LIBMESH_DIM) i = 0; i < LIBMESH_DIM; i++)
226 14916 : translated_node(i) = input_node(i);
227 44 : }
228 : }
229 :
230 42 : final_mesh->set_isnt_prepared();
231 42 : return dynamic_pointer_cast<MeshBase>(final_mesh);
232 42 : }
233 : }
234 :
235 : void
236 2134 : CombinerGenerator::copyIntoMesh(UnstructuredMesh & destination, const UnstructuredMesh & source)
237 : {
238 2134 : dof_id_type node_delta = destination.max_node_id();
239 2134 : dof_id_type elem_delta = destination.max_elem_id();
240 :
241 : unique_id_type unique_delta =
242 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
243 2134 : destination.parallel_max_unique_id();
244 : #else
245 : 0;
246 : #endif
247 :
248 : // Prevent overlaps by offsetting the subdomains in
249 2134 : std::unordered_map<subdomain_id_type, subdomain_id_type> id_remapping;
250 2134 : unsigned int block_offset = 0;
251 2134 : if (_avoid_merging_subdomains)
252 : {
253 : // Note: if performance becomes an issue, this is overkill for just getting the max node id
254 147 : std::set<subdomain_id_type> source_ids;
255 147 : std::set<subdomain_id_type> dest_ids;
256 147 : source.subdomain_ids(source_ids, true);
257 147 : destination.subdomain_ids(dest_ids, true);
258 : mooseAssert(source_ids.size(), "Should have a subdomain");
259 : mooseAssert(dest_ids.size(), "Should have a subdomain");
260 147 : unsigned int max_dest_bid = *dest_ids.rbegin();
261 147 : unsigned int min_source_bid = *source_ids.begin();
262 147 : _communicator.max(max_dest_bid);
263 147 : _communicator.min(min_source_bid);
264 147 : block_offset = 1 + max_dest_bid - min_source_bid;
265 294 : for (const auto bid : source_ids)
266 147 : id_remapping[bid] = block_offset + bid;
267 147 : }
268 :
269 : // Copy mesh data over from the other mesh
270 2134 : destination.copy_nodes_and_elements(source,
271 : // Skipping this should cause the neighbors
272 : // to simply be copied from the other mesh
273 : // (which makes sense and is way faster)
274 : /*skip_find_neighbors = */ true,
275 : elem_delta,
276 : node_delta,
277 : unique_delta,
278 2134 : _avoid_merging_subdomains ? &id_remapping : nullptr);
279 :
280 : // Get an offset to prevent overlaps / wild merging between boundaries
281 2134 : BoundaryInfo & boundary = destination.get_boundary_info();
282 2134 : const BoundaryInfo & other_boundary = source.get_boundary_info();
283 :
284 2134 : unsigned int bid_offset = 0;
285 2134 : if (_avoid_merging_boundaries)
286 : {
287 180 : const auto boundary_ids = boundary.get_boundary_ids();
288 180 : const auto other_boundary_ids = other_boundary.get_boundary_ids();
289 180 : unsigned int max_dest_bid = boundary_ids.size() ? *boundary_ids.rbegin() : 0;
290 180 : unsigned int min_source_bid = other_boundary_ids.size() ? *other_boundary_ids.begin() : 0;
291 180 : _communicator.max(max_dest_bid);
292 180 : _communicator.min(min_source_bid);
293 180 : bid_offset = 1 + max_dest_bid - min_source_bid;
294 180 : }
295 :
296 : // Note: the code below originally came from ReplicatedMesh::stitch_mesh_helper()
297 : // in libMesh replicated_mesh.C around line 1203
298 :
299 : // Copy BoundaryInfo from other_mesh too. We do this via the
300 : // list APIs rather than element-by-element for speed.
301 60172 : for (const auto & t : other_boundary.build_node_list())
302 60172 : boundary.add_node(std::get<0>(t) + node_delta, bid_offset + std::get<1>(t));
303 :
304 40845 : for (const auto & t : other_boundary.build_side_list())
305 40845 : boundary.add_side(std::get<0>(t) + elem_delta, std::get<1>(t), bid_offset + std::get<2>(t));
306 :
307 2134 : for (const auto & t : other_boundary.build_edge_list())
308 2134 : boundary.add_edge(std::get<0>(t) + elem_delta, std::get<1>(t), bid_offset + std::get<2>(t));
309 :
310 2134 : for (const auto & t : other_boundary.build_shellface_list())
311 0 : boundary.add_shellface(
312 2134 : std::get<0>(t) + elem_delta, std::get<1>(t), bid_offset + std::get<2>(t));
313 :
314 : // Check for the case with two block ids sharing the same name
315 2134 : if (_avoid_merging_subdomains)
316 : {
317 327 : for (const auto & [block_id, block_name] : destination.get_subdomain_name_map())
318 236 : for (const auto & [source_id, source_name] : source.get_subdomain_name_map())
319 56 : if (block_name == source_name)
320 0 : paramWarning("avoid_merging_subdomains",
321 0 : "Not merging subdomains is creating two subdomains with the same name '" +
322 0 : block_name + "' but different ids: " + std::to_string(source_id) +
323 0 : " & " + std::to_string(block_id + block_offset) +
324 : ".\n We recommend using a RenameBlockGenerator to prevent this as you "
325 : "will get errors reading the Exodus output later.");
326 : }
327 :
328 2291 : for (const auto & [block_id, block_name] : source.get_subdomain_name_map())
329 314 : destination.set_subdomain_name_map().insert(
330 314 : std::make_pair<SubdomainID, SubdomainName>(block_id + block_offset, block_name));
331 :
332 : // Check for the case with two boundary ids sharing the same name
333 2134 : if (_avoid_merging_boundaries)
334 : {
335 922 : for (const auto & [b_id, b_name] : other_boundary.get_sideset_name_map())
336 5266 : for (const auto & [source_id, source_name] : boundary.get_sideset_name_map())
337 4524 : if (b_name == source_name)
338 0 : paramWarning(
339 : "avoid_merging_boundaries",
340 0 : "Not merging boundaries is creating two sidesets with the same name '" + b_name +
341 0 : "' but different ids: " + std::to_string(source_id) + " & " +
342 0 : std::to_string(b_id + bid_offset) +
343 : ".\n We recommend using a RenameBoundaryGenerator to prevent this as you "
344 : "will get errors reading the Exodus output later.");
345 922 : for (const auto & [b_id, b_name] : other_boundary.get_nodeset_name_map())
346 5266 : for (const auto & [source_id, source_name] : boundary.get_nodeset_name_map())
347 4524 : if (b_name == source_name)
348 0 : paramWarning(
349 : "avoid_merging_boundaries",
350 0 : "Not merging boundaries is creating two nodesets with the same name '" + b_name +
351 0 : "' but different ids: " + std::to_string(source_id) + " & " +
352 0 : std::to_string(b_id + bid_offset) +
353 : ".\n We recommend using a RenameBoundaryGenerator to prevent this as you "
354 : "will get errors reading the Exodus output later.");
355 : }
356 :
357 9760 : for (const auto & [nodeset_id, nodeset_name] : other_boundary.get_nodeset_name_map())
358 15252 : boundary.set_nodeset_name_map().insert(
359 15252 : std::make_pair<BoundaryID, BoundaryName>(nodeset_id + bid_offset, nodeset_name));
360 :
361 10046 : for (const auto & [sideset_id, sideset_name] : other_boundary.get_sideset_name_map())
362 15824 : boundary.set_sideset_name_map().insert(
363 15824 : std::make_pair<BoundaryID, BoundaryName>(sideset_id + bid_offset, sideset_name));
364 :
365 2134 : for (const auto & [edgeset_id, edgeset_name] : other_boundary.get_edgeset_name_map())
366 0 : boundary.set_edgeset_name_map().insert(
367 0 : std::make_pair<BoundaryID, BoundaryName>(edgeset_id + bid_offset, edgeset_name));
368 2134 : }
|
