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