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 15501 : CombinerGenerator::validParams()
26 : {
27 15501 : InputParameters params = MeshGenerator::validParams();
28 :
29 15501 : 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 15501 : 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 15501 : 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 15501 : params.addParam<std::vector<FileName>>(
45 : "positions_file", "Alternative way to provide the position of each given mesh.");
46 :
47 46503 : params.addParam<bool>("avoid_merging_subdomains",
48 31002 : 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 46503 : params.addParam<bool>("avoid_merging_boundaries",
53 31002 : 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 15501 : return params;
59 0 : }
60 :
61 618 : CombinerGenerator::CombinerGenerator(const InputParameters & parameters)
62 : : MeshGenerator(parameters),
63 618 : _meshes(getMeshes("inputs")),
64 618 : _input_names(getParam<std::vector<MeshGeneratorName>>("inputs")),
65 618 : _avoid_merging_subdomains(getParam<bool>("avoid_merging_subdomains")),
66 1854 : _avoid_merging_boundaries(getParam<bool>("avoid_merging_boundaries"))
67 : {
68 618 : if (_input_names.empty())
69 0 : paramError("input_names", "You need to specify at least one MeshGenerator as an input.");
70 :
71 618 : 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 618 : if (_input_names.size() == 1)
77 48 : 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 618 : }
82 :
83 : void
84 601 : CombinerGenerator::fillPositions()
85 : {
86 601 : if (isParamValid("positions"))
87 : {
88 241 : _positions = getParam<std::vector<Point>>("positions");
89 :
90 : // the check in the constructor wont catch error where the user sets positions = ''
91 241 : 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 237 : if (_input_names.size() != 1)
95 208 : 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 360 : else if (isParamValid("positions_file"))
102 : {
103 28 : 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 28 : 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 44 : for (const auto & f : positions_file)
111 : {
112 24 : MooseUtils::DelimitedFileReader file(f, &_communicator);
113 24 : file.setFormatFlag(MooseUtils::DelimitedFileReader::FormatFlag::ROWS);
114 24 : file.read();
115 :
116 24 : const std::vector<Point> & data = file.getDataAsPoints();
117 :
118 24 : if (_input_names.size() != 1)
119 14 : 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 80 : for (const auto & d : data)
125 60 : _positions.push_back(d);
126 20 : }
127 20 : }
128 585 : }
129 :
130 : std::unique_ptr<MeshBase>
131 601 : CombinerGenerator::generate()
132 : {
133 : // Two cases:
134 : // 1. Multiple input meshes and optional positions
135 : // 2. One input mesh and multiple positions
136 601 : fillPositions();
137 :
138 : // Case 1
139 585 : if (_meshes.size() != 1)
140 : {
141 : // merge all meshes into the first one
142 546 : auto mesh = dynamic_pointer_cast<UnstructuredMesh>(*_meshes[0]);
143 :
144 546 : 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 546 : if (_positions.size())
149 : {
150 214 : MeshTools::Modification::translate(
151 214 : *mesh, _positions[0](0), _positions[0](1), _positions[0](2));
152 : }
153 :
154 : // Read in all of the other meshes
155 2420 : for (MooseIndex(_meshes) i = 1; i < _meshes.size(); ++i)
156 : {
157 1874 : auto other_mesh = dynamic_pointer_cast<UnstructuredMesh>(*_meshes[i]);
158 :
159 1874 : if (!other_mesh)
160 0 : paramError("inputs", _input_names[i], " is not a valid unstructured mesh");
161 :
162 : // Move It
163 1874 : if (_positions.size())
164 : {
165 698 : MeshTools::Modification::translate(
166 698 : *other_mesh, _positions[i](0), _positions[i](1), _positions[i](2));
167 : }
168 :
169 1874 : copyIntoMesh(*mesh, *other_mesh);
170 1874 : }
171 :
172 546 : mesh->set_isnt_prepared();
173 546 : return dynamic_pointer_cast<MeshBase>(mesh);
174 546 : }
175 : else // Case 2
176 : {
177 39 : auto input_mesh = dynamic_pointer_cast<UnstructuredMesh>(*_meshes[0]);
178 :
179 39 : 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 39 : input_mesh->clone(); // This is required because dynamic_pointer_cast() requires an l-value
185 39 : auto final_mesh = dynamic_pointer_cast<UnstructuredMesh>(copy);
186 :
187 39 : if (!final_mesh)
188 0 : mooseError("Unable to copy mesh!");
189 :
190 39 : MeshTools::Modification::translate(
191 39 : *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 39 : copy = input_mesh->clone();
205 39 : auto translated_mesh = dynamic_pointer_cast<UnstructuredMesh>(copy);
206 :
207 39 : if (!translated_mesh)
208 0 : mooseError("Unable to copy mesh!");
209 :
210 79 : for (MooseIndex(_meshes) i = 1; i < _positions.size(); ++i)
211 : {
212 : // Move
213 40 : MeshTools::Modification::translate(
214 40 : *translated_mesh, _positions[i](0), _positions[i](1), _positions[i](2));
215 :
216 : // Copy into final mesh
217 40 : copyIntoMesh(*final_mesh, *translated_mesh);
218 :
219 : // Reset nodal coordinates
220 9016 : for (auto translated_node_ptr : translated_mesh->node_ptr_range())
221 : {
222 4488 : auto & translated_node = *translated_node_ptr;
223 4488 : auto & input_node = input_mesh->node_ref(translated_node_ptr->id());
224 :
225 17952 : for (MooseIndex(LIBMESH_DIM) i = 0; i < LIBMESH_DIM; i++)
226 13464 : translated_node(i) = input_node(i);
227 40 : }
228 : }
229 :
230 39 : final_mesh->set_isnt_prepared();
231 39 : return dynamic_pointer_cast<MeshBase>(final_mesh);
232 39 : }
233 : }
234 :
235 : void
236 1914 : CombinerGenerator::copyIntoMesh(UnstructuredMesh & destination, const UnstructuredMesh & source)
237 : {
238 1914 : dof_id_type node_delta = destination.max_node_id();
239 1914 : dof_id_type elem_delta = destination.max_elem_id();
240 :
241 : unique_id_type unique_delta =
242 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
243 1914 : destination.parallel_max_unique_id();
244 : #else
245 : 0;
246 : #endif
247 :
248 : // Prevent overlaps by offsetting the subdomains in
249 1914 : std::unordered_map<subdomain_id_type, subdomain_id_type> id_remapping;
250 1914 : unsigned int block_offset = 0;
251 1914 : if (_avoid_merging_subdomains)
252 : {
253 : // Note: if performance becomes an issue, this is overkill for just getting the max node id
254 136 : std::set<subdomain_id_type> source_ids;
255 136 : std::set<subdomain_id_type> dest_ids;
256 136 : source.subdomain_ids(source_ids, true);
257 136 : 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 136 : unsigned int max_dest_bid = *dest_ids.rbegin();
261 136 : unsigned int min_source_bid = *source_ids.begin();
262 136 : _communicator.max(max_dest_bid);
263 136 : _communicator.min(min_source_bid);
264 136 : block_offset = 1 + max_dest_bid - min_source_bid;
265 272 : for (const auto bid : source_ids)
266 136 : id_remapping[bid] = block_offset + bid;
267 136 : }
268 :
269 : // Copy mesh data over from the other mesh
270 1914 : 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 1914 : _avoid_merging_subdomains ? &id_remapping : nullptr);
279 :
280 : // Get an offset to prevent overlaps / wild merging between boundaries
281 1914 : BoundaryInfo & boundary = destination.get_boundary_info();
282 1914 : const BoundaryInfo & other_boundary = source.get_boundary_info();
283 :
284 1914 : unsigned int bid_offset = 0;
285 1914 : if (_avoid_merging_boundaries)
286 : {
287 166 : const auto boundary_ids = boundary.get_boundary_ids();
288 166 : const auto other_boundary_ids = other_boundary.get_boundary_ids();
289 166 : unsigned int max_dest_bid = boundary_ids.size() ? *boundary_ids.rbegin() : 0;
290 166 : unsigned int min_source_bid = other_boundary_ids.size() ? *other_boundary_ids.begin() : 0;
291 166 : _communicator.max(max_dest_bid);
292 166 : _communicator.min(min_source_bid);
293 166 : bid_offset = 1 + max_dest_bid - min_source_bid;
294 166 : }
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 54638 : for (const auto & t : other_boundary.build_node_list())
302 54638 : boundary.add_node(std::get<0>(t) + node_delta, bid_offset + std::get<1>(t));
303 :
304 36956 : for (const auto & t : other_boundary.build_side_list())
305 36956 : boundary.add_side(std::get<0>(t) + elem_delta, std::get<1>(t), bid_offset + std::get<2>(t));
306 :
307 1914 : for (const auto & t : other_boundary.build_edge_list())
308 1914 : boundary.add_edge(std::get<0>(t) + elem_delta, std::get<1>(t), bid_offset + std::get<2>(t));
309 :
310 1914 : for (const auto & t : other_boundary.build_shellface_list())
311 0 : boundary.add_shellface(
312 1914 : 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 1914 : if (_avoid_merging_subdomains)
316 : {
317 303 : for (const auto & [block_id, block_name] : destination.get_subdomain_name_map())
318 219 : for (const auto & [source_id, source_name] : source.get_subdomain_name_map())
319 52 : 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 2029 : for (const auto & [block_id, block_name] : source.get_subdomain_name_map())
329 230 : destination.set_subdomain_name_map().insert(
330 230 : 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 1914 : if (_avoid_merging_boundaries)
334 : {
335 850 : for (const auto & [b_id, b_name] : other_boundary.get_sideset_name_map())
336 4852 : for (const auto & [source_id, source_name] : boundary.get_sideset_name_map())
337 4168 : 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 850 : for (const auto & [b_id, b_name] : other_boundary.get_nodeset_name_map())
346 4852 : for (const auto & [source_id, source_name] : boundary.get_nodeset_name_map())
347 4168 : 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 8824 : for (const auto & [nodeset_id, nodeset_name] : other_boundary.get_nodeset_name_map())
358 13820 : boundary.set_nodeset_name_map().insert(
359 13820 : std::make_pair<BoundaryID, BoundaryName>(nodeset_id + bid_offset, nodeset_name));
360 :
361 9084 : for (const auto & [sideset_id, sideset_name] : other_boundary.get_sideset_name_map())
362 14340 : boundary.set_sideset_name_map().insert(
363 14340 : std::make_pair<BoundaryID, BoundaryName>(sideset_id + bid_offset, sideset_name));
364 :
365 1914 : 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 1914 : }
|
