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 "PatchSidesetGenerator.h"
11 : #include "InputParameters.h"
12 : #include "MooseTypes.h"
13 : #include "CastUniquePointer.h"
14 : #include "MooseUtils.h"
15 : #include "MooseMeshUtils.h"
16 :
17 : #include "libmesh/distributed_mesh.h"
18 : #include "libmesh/elem.h"
19 : #include "libmesh/linear_partitioner.h"
20 : #include "libmesh/centroid_partitioner.h"
21 : #include "libmesh/parmetis_partitioner.h"
22 : #include "libmesh/hilbert_sfc_partitioner.h"
23 : #include "libmesh/morton_sfc_partitioner.h"
24 : #include "libmesh/enum_elem_type.h"
25 :
26 : #include <set>
27 : #include <limits>
28 : #include "libmesh/mesh_tools.h"
29 :
30 : registerMooseObject("HeatTransferApp", PatchSidesetGenerator);
31 :
32 : InputParameters
33 1054 : PatchSidesetGenerator::validParams()
34 : {
35 1054 : InputParameters params = MeshGenerator::validParams();
36 :
37 2108 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
38 2108 : params.addRequiredParam<BoundaryName>("boundary",
39 : "The boundary that will be divided into patches");
40 2108 : params.addRequiredRangeCheckedParam<unsigned int>(
41 : "n_patches", "n_patches>0", "Number of patches");
42 :
43 1054 : MooseEnum partitioning = MooseMesh::partitioning(); // default MOOSE partitioning
44 1054 : partitioning += "grid"; // ...but also add our own
45 2108 : params.addParam<MooseEnum>(
46 : "partitioner",
47 : partitioning,
48 : "Specifies a mesh partitioner to use when splitting the mesh for a parallel computation.");
49 :
50 2108 : MooseEnum direction("x y z radial");
51 2108 : params.addParam<MooseEnum>("centroid_partitioner_direction",
52 : direction,
53 : "Specifies the sort direction if using the centroid partitioner. "
54 : "Available options: x, y, z, radial");
55 :
56 2108 : params.addParamNamesToGroup("partitioner centroid_partitioner_direction", "Partitioning");
57 :
58 1054 : params.addClassDescription(
59 : "Divides the given sideset into smaller patches of roughly equal size.");
60 :
61 1054 : return params;
62 1054 : }
63 :
64 527 : PatchSidesetGenerator::PatchSidesetGenerator(const InputParameters & parameters)
65 : : MeshGenerator(parameters),
66 527 : _input(getMesh("input")),
67 1054 : _n_patches(getParam<unsigned int>("n_patches")),
68 1054 : _sideset_name(getParam<BoundaryName>("boundary")),
69 1581 : _partitioner_name(getParam<MooseEnum>("partitioner"))
70 : {
71 527 : }
72 :
73 : std::unique_ptr<MeshBase>
74 392 : PatchSidesetGenerator::generate()
75 : {
76 392 : std::unique_ptr<MeshBase> mesh = std::move(_input);
77 :
78 784 : _mesh->errorIfDistributedMesh("PatchSidesetGenerator");
79 :
80 : // Get a reference to our BoundaryInfo object for later use
81 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
82 :
83 : // get dimensionality
84 392 : _dim = mesh->mesh_dimension() - 1;
85 :
86 : // get a list of all sides; vector of tuples (elem, loc_side, side_set)
87 392 : auto side_list = boundary_info.build_active_side_list();
88 :
89 : // check if the provided sideset name is actually a sideset id
90 : // if _sideset_name can be converted to integer it's interpreted
91 : // as sideset id
92 392 : std::stringstream ss;
93 392 : ss << _sideset_name;
94 392 : ss >> _sideset;
95 392 : if (ss.fail())
96 88 : _sideset = boundary_info.get_id_by_name(_sideset_name);
97 :
98 : // make sure that _sideset exists
99 392 : if (_sideset == BoundaryInfo::invalid_id)
100 0 : paramError("sideset_name", "Not a valid boundary");
101 :
102 : // create a dim - 1 dimensional mesh
103 392 : auto boundary_mesh = buildReplicatedMesh(mesh->mesh_dimension() - 1);
104 784 : boundary_mesh->set_mesh_dimension(mesh->mesh_dimension() - 1);
105 392 : boundary_mesh->set_spatial_dimension(mesh->mesh_dimension());
106 :
107 : // nodes in the new mesh by boundary_node_id (index)
108 : std::vector<Node *> boundary_nodes;
109 : // a map from the node numbering on the volumetric mesh to the numbering
110 : // on the boundary_mesh
111 : std::map<dof_id_type, dof_id_type> mesh_node_id_to_boundary_node_id;
112 : // a local counter keeping track of how many entries have been added to boundary_nodes
113 : dof_id_type boundary_node_id = 0;
114 : // a map from new element id in the boundary mesh to the element id/side/sideset
115 : // tuple it came from
116 : std::map<dof_id_type, std::tuple<dof_id_type, unsigned short int, boundary_id_type>>
117 : boundary_elem_to_mesh_elem;
118 115820 : for (auto & side : side_list)
119 : {
120 115428 : if (std::get<2>(side) == _sideset)
121 : {
122 : // the original volumetric mesh element
123 16206 : const Elem * elem = mesh->elem_ptr(std::get<0>(side));
124 :
125 : // the boundary element
126 16206 : std::unique_ptr<const Elem> boundary_elem = elem->side_ptr(std::get<1>(side));
127 :
128 : // an array that saves the boundary node ids of this elem in the right order
129 16206 : std::vector<dof_id_type> bnd_elem_node_ids(boundary_elem->n_nodes());
130 :
131 : // loop through the nodes in boundary_elem
132 78038 : for (MooseIndex(boundary_elem->n_nodes()) j = 0; j < boundary_elem->n_nodes(); ++j)
133 : {
134 : const Node * node = boundary_elem->node_ptr(j);
135 :
136 : // Is this node a new node?
137 61832 : if (mesh_node_id_to_boundary_node_id.find(node->id()) ==
138 : mesh_node_id_to_boundary_node_id.end())
139 : {
140 : // yes, it is new, need to add it to the mesh_node_id_to_boundary_node_id map
141 18346 : mesh_node_id_to_boundary_node_id.insert(
142 18346 : std::pair<dof_id_type, dof_id_type>(node->id(), boundary_node_id));
143 :
144 : // this adds this node to the boundary mesh and puts it at the right position
145 : // in the boundary_nodes array
146 18346 : Point pt(*node);
147 18346 : boundary_nodes.push_back(boundary_mesh->add_point(pt, boundary_node_id));
148 :
149 : // keep track of the boundary node for setting up the element
150 18346 : bnd_elem_node_ids[j] = boundary_node_id;
151 :
152 : // increment the boundary_node_id counter
153 18346 : ++boundary_node_id;
154 : }
155 : else
156 43486 : bnd_elem_node_ids[j] = mesh_node_id_to_boundary_node_id.find(node->id())->second;
157 : }
158 :
159 : // all nodes for this element have been added, so we can add the element to the
160 : // boundary mesh
161 16206 : Elem * new_bnd_elem = boundaryElementHelper(*boundary_mesh, boundary_elem->type());
162 :
163 : // keep track of these new boundary elements in boundary_elem_to_mesh_elem
164 16206 : boundary_elem_to_mesh_elem.insert(
165 16206 : std::pair<dof_id_type, std::tuple<dof_id_type, unsigned short int, boundary_id_type>>(
166 : new_bnd_elem->id(), side));
167 :
168 : // set the nodes & subdomain_id of the new element by looping over the
169 : // boundary_elem and then inserting its nodes into new_bnd_elem in the
170 : // same order
171 78038 : for (MooseIndex(boundary_elem->n_nodes()) j = 0; j < boundary_elem->n_nodes(); ++j)
172 : {
173 61832 : dof_id_type old_node_id = boundary_elem->node_ptr(j)->id();
174 61832 : if (mesh_node_id_to_boundary_node_id.find(old_node_id) ==
175 : mesh_node_id_to_boundary_node_id.end())
176 0 : mooseError("Node id", old_node_id, " not linked to new node id.");
177 61832 : dof_id_type new_node_id = mesh_node_id_to_boundary_node_id.find(old_node_id)->second;
178 61832 : new_bnd_elem->set_node(j, boundary_nodes[new_node_id]);
179 : }
180 16206 : }
181 : }
182 :
183 : // partition the boundary mesh
184 392 : boundary_mesh->prepare_for_use();
185 392 : _n_boundary_mesh_elems = boundary_mesh->n_elem();
186 392 : if (_partitioner_name == "grid")
187 15 : partition(*boundary_mesh);
188 : else
189 : {
190 754 : auto partitioner_enum = getParam<MooseEnum>("partitioner");
191 377 : MooseMesh::setPartitioner(*boundary_mesh, partitioner_enum, false, _pars, *this);
192 375 : boundary_mesh->partition(_n_patches);
193 375 : }
194 :
195 : // make sure every partition has at least one element; if not rename and adjust _n_patches
196 390 : checkPartitionAndCompress(*boundary_mesh);
197 :
198 : // prepare sideset names and boundary_ids added to mesh
199 : std::vector<BoundaryName> sideset_names =
200 390 : sidesetNameHelper(boundary_info.get_sideset_name(_sideset));
201 :
202 : std::vector<boundary_id_type> boundary_ids =
203 390 : MooseMeshUtils::getBoundaryIDs(*mesh, sideset_names, true);
204 :
205 : mooseAssert(sideset_names.size() == _n_patches,
206 : "sideset_names must have as many entries as user-requested number of patches.");
207 : mooseAssert(boundary_ids.size() == _n_patches,
208 : "boundary_ids must have as many entries as user-requested number of patches.");
209 :
210 : // loop through all elements in the boundary mesh and assign the side of
211 : // the _original_ element to the new sideset
212 30692 : for (const auto & elem : boundary_mesh->active_element_ptr_range())
213 : {
214 14956 : if (boundary_elem_to_mesh_elem.find(elem->id()) == boundary_elem_to_mesh_elem.end())
215 0 : mooseError("Element in the boundary mesh with id ",
216 0 : elem->id(),
217 : " not found in boundary_elem_to_mesh_elem.");
218 :
219 14956 : auto side = boundary_elem_to_mesh_elem.find(elem->id())->second;
220 :
221 : mooseAssert(elem->processor_id() < boundary_ids.size(),
222 : "Processor id larger than number of patches.");
223 14956 : boundary_info.add_side(
224 14956 : std::get<0>(side), std::get<1>(side), boundary_ids[elem->processor_id()]);
225 390 : }
226 :
227 : // make sure new boundary names are set
228 1519 : for (MooseIndex(boundary_ids.size()) j = 0; j < boundary_ids.size(); ++j)
229 : {
230 1129 : boundary_info.sideset_name(boundary_ids[j]) = sideset_names[j];
231 1129 : boundary_info.nodeset_name(boundary_ids[j]) = sideset_names[j];
232 : }
233 :
234 390 : return mesh;
235 780 : }
236 :
237 : void
238 15 : PatchSidesetGenerator::partition(MeshBase & mesh)
239 : {
240 15 : if (_partitioner_name == "grid")
241 : {
242 : // Figure out the physical bounds of the given mesh
243 15 : auto bounding_box = MeshTools::create_bounding_box(mesh);
244 : const auto & min = bounding_box.min();
245 : const auto & max = bounding_box.max();
246 : const auto & delta = max - min;
247 :
248 : // set number of elements
249 15 : std::vector<unsigned int> nelems(3);
250 15 : if (_dim == 1)
251 : {
252 : // find the largest component in delta
253 : unsigned int largest_id = 0;
254 : Real largest = delta(0);
255 30 : for (unsigned int j = 1; j < 3; ++j)
256 20 : if (largest < delta(j))
257 : {
258 : largest = delta(j);
259 : largest_id = j;
260 : }
261 :
262 : // set nelems now
263 10 : nelems = {1, 1, 1};
264 10 : nelems[largest_id] = _n_patches;
265 : }
266 : else
267 : {
268 : // find the smallest component in delta
269 : unsigned int smallest_id = 0;
270 : Real smallest = delta(0);
271 15 : for (unsigned int j = 1; j < 3; ++j)
272 10 : if (smallest > delta(j))
273 : {
274 : smallest = delta(j);
275 : smallest_id = j;
276 : }
277 :
278 : // store the ids for the two larger dimensions
279 : unsigned int id1 = 1, id2 = 2;
280 5 : if (smallest_id == 1)
281 : id1 = 0;
282 5 : else if (smallest_id == 2)
283 : id2 = 0;
284 :
285 : // set number of elements
286 5 : nelems[smallest_id] = 1;
287 5 : nelems[id1] = std::round(std::sqrt(delta(id1) / delta(id2) * _n_patches));
288 5 : nelems[id2] = std::round(std::sqrt(delta(id2) / delta(id1) * _n_patches));
289 5 : unsigned int final_n_patches = nelems[id1] * nelems[id2];
290 : // We need to check if the number of requested patches and the number of
291 : // actually created patches matches. If the two do not match, then a warning
292 : // is printed.
293 5 : if (_n_patches != final_n_patches)
294 : {
295 5 : _console << "Note: For creating radiation patches for boundary " << _sideset
296 5 : << " using grid partitioner number of patches was changed from " << _n_patches
297 5 : << " to " << final_n_patches << std::endl;
298 5 : _n_patches = final_n_patches;
299 : }
300 : }
301 :
302 15 : const Real dx = delta(0) / nelems[0];
303 15 : const Real dy = delta(1) / nelems[1];
304 15 : const Real dz = delta(2) / nelems[2];
305 3755 : for (auto & elem_ptr : mesh.active_element_ptr_range())
306 : {
307 3725 : const Point centroid = elem_ptr->vertex_average();
308 : processor_id_type proc_id;
309 3725 : const unsigned int ix = std::floor((centroid(0) - min(0)) / dx);
310 3725 : const unsigned int iy = std::floor((centroid(1) - min(1)) / dy);
311 3725 : const unsigned int iz = std::floor((centroid(2) - min(2)) / dz);
312 3725 : proc_id = ix + iy * nelems[0] + iz * nelems[0] * nelems[1];
313 3725 : elem_ptr->processor_id() = proc_id;
314 15 : }
315 15 : }
316 : else
317 0 : mooseError("Partitioner ", _partitioner_name, " not recognized.");
318 15 : }
319 :
320 : void
321 390 : PatchSidesetGenerator::checkPartitionAndCompress(MeshBase & mesh)
322 : {
323 : std::set<processor_id_type> processor_ids;
324 30692 : for (auto & elem_ptr : mesh.active_element_ptr_range())
325 15346 : processor_ids.insert(elem_ptr->processor_id());
326 :
327 390 : if (processor_ids.size() == _n_patches)
328 : return;
329 :
330 : // at least one partition does not have an elem assigned to it
331 : // adjust _n_patches
332 5 : _console << "Some partitions for side set " << _sideset
333 5 : << " are empty. Adjusting number of patches from " << _n_patches << " to "
334 5 : << processor_ids.size() << std::endl;
335 5 : _n_patches = processor_ids.size();
336 :
337 : // create a vector and sort it
338 : std::vector<processor_id_type> processor_ids_vec;
339 255 : for (auto & p : processor_ids)
340 250 : processor_ids_vec.push_back(p);
341 5 : std::sort(processor_ids_vec.begin(), processor_ids_vec.end());
342 :
343 : // now remap the processor ids
344 : std::map<processor_id_type, processor_id_type> processor_id_remap;
345 255 : for (MooseIndex(processor_ids_vec.size()) j = 0; j < processor_ids_vec.size(); ++j)
346 250 : processor_id_remap[processor_ids_vec[j]] = j;
347 :
348 310 : for (auto & elem_ptr : mesh.active_element_ptr_range())
349 : {
350 300 : processor_id_type p = elem_ptr->processor_id();
351 : const auto & it = processor_id_remap.find(p);
352 300 : if (it == processor_id_remap.end())
353 0 : mooseError("Parition id ", p, " not in processor_id_remap.");
354 300 : elem_ptr->processor_id() = it->second;
355 5 : }
356 5 : }
357 :
358 : std::vector<BoundaryName>
359 390 : PatchSidesetGenerator::sidesetNameHelper(const std::string & base_name) const
360 : {
361 : std::vector<BoundaryName> rv;
362 1519 : for (unsigned int j = 0; j < _n_patches; ++j)
363 : {
364 1129 : std::stringstream ss;
365 1129 : ss << base_name << "_" << j;
366 1129 : rv.push_back(ss.str());
367 1129 : }
368 390 : return rv;
369 0 : }
370 :
371 : Elem *
372 16206 : PatchSidesetGenerator::boundaryElementHelper(MeshBase & mesh, libMesh::ElemType type) const
373 : {
374 16206 : std::unique_ptr<Elem> elem = libMesh::Elem::build(type);
375 16206 : if (elem->dim() < 3)
376 32412 : return mesh.add_elem(std::move(elem));
377 :
378 0 : mooseError("Unsupported element type (libMesh elem_type enum): ", type);
379 16206 : }
|
