Line data Source code
1 : /********************************************************************/
2 : /* SOFTWARE COPYRIGHT NOTIFICATION */
3 : /* Cardinal */
4 : /* */
5 : /* (c) 2021 UChicago Argonne, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /* */
8 : /* Prepared by UChicago Argonne, LLC */
9 : /* Under Contract No. DE-AC02-06CH11357 */
10 : /* With the U. S. Department of Energy */
11 : /* */
12 : /* Prepared by Battelle Energy Alliance, LLC */
13 : /* Under Contract No. DE-AC07-05ID14517 */
14 : /* With the U. S. Department of Energy */
15 : /* */
16 : /* See LICENSE for full restrictions */
17 : /********************************************************************/
18 :
19 : #include "NekMeshGenerator.h"
20 : #include "CastUniquePointer.h"
21 : #include "UserErrorChecking.h"
22 : #include "MooseMeshUtils.h"
23 : #include "DelimitedFileReader.h"
24 : #include "GeometryUtils.h"
25 :
26 : #include "libmesh/mesh_tools.h"
27 : #include "libmesh/cell_hex8.h"
28 : #include "libmesh/cell_hex20.h"
29 : #include "libmesh/cell_hex27.h"
30 : #include "libmesh/face_quad4.h"
31 : #include "libmesh/face_quad8.h"
32 : #include "libmesh/face_quad9.h"
33 :
34 : registerMooseObject("CardinalApp", NekMeshGenerator);
35 :
36 : InputParameters
37 894 : NekMeshGenerator::validParams()
38 : {
39 894 : InputParameters params = MeshGenerator::validParams();
40 1788 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
41 :
42 1788 : MooseEnum geom("cylinder sphere");
43 1788 : params.addParam<MooseEnum>(
44 : "geometry_type", geom, "Geometry type to use for moving boundary nodes");
45 :
46 : // optional parameters if fitting sidesets to a circular surface
47 1788 : MooseEnum axis("x y z", "z");
48 1788 : params.addParam<MooseEnum>(
49 : "axis",
50 : axis,
51 : "If 'geometry_type = cylinder', the axis of the mesh about which to build "
52 : "the cylinder surface(s)");
53 1788 : params.addParam<std::vector<BoundaryName>>("boundary",
54 : "Boundary(s) to enforce the curved surface");
55 1788 : params.addParam<std::vector<Real>>("radius", "Radius(es) of the surfaces");
56 1788 : params.addParam<std::vector<std::vector<Real>>>(
57 : "origins",
58 : "Origin(s) about which to form the curved surfaces; "
59 : "if not specified, all values default to (0, 0, 0)");
60 1788 : params.addParam<std::vector<std::string>>("origins_files",
61 : "Origin(s) about which to form the curved surfaces, "
62 : "with a file of points provided for each boundary. If "
63 : "not specified, all values default to (0, 0, 0)");
64 1788 : params.addParam<std::vector<unsigned int>>(
65 : "layers",
66 : "Number of layers to sweep for each "
67 : "boundary when forming the curved surfaces; if not specified, all values default to 0");
68 :
69 : // optional parameters if fitting corners of a polygon to radius of curvature
70 1788 : params.addParam<bool>("curve_corners",
71 1788 : false,
72 : "If 'geometry_type = cylinder', whether to move elements to respect radius "
73 : "of curvature of polygon corners");
74 1788 : params.addRangeCheckedParam<unsigned int>(
75 : "polygon_sides",
76 : "polygon_sides > 2",
77 : "If 'geometry_type = cylinder' and when curving corners, the number of sides of the polygon "
78 : "to use for identifying corners");
79 1788 : params.addRangeCheckedParam<Real>("polygon_size",
80 : "polygon_size > 0.0",
81 : "If 'geometry_type = cylinder' and when curving corners, the "
82 : "size of the polygon (measured as distance from center to a "
83 : "corner) to use for identifying corners");
84 1788 : params.addRangeCheckedParam<Real>("corner_radius",
85 : "corner_radius >= 0.0",
86 : "If 'geometry_type = cylinder' and when curving corners, the "
87 : "radius of curvature of the corners");
88 1788 : params.addParam<unsigned int>("polygon_layers",
89 1788 : 0,
90 : "If 'geometry_type = cylinder' and when curving corners, the "
91 : "number of layers to sweep for each polygon corner");
92 1788 : params.addParam<std::vector<Real>>(
93 : "polygon_layer_smoothing",
94 : "If 'geometry_type = cylinder' and when curving corners, the multiplicative factor to apply "
95 : "to each boundary layer; if not "
96 : "specified, all values default to 1.0");
97 1788 : params.addParam<BoundaryName>(
98 : "polygon_boundary",
99 : "If 'geometry_type = cylinder', boundary to enforce radius of curvature for polygon corners");
100 1788 : params.addParam<std::vector<std::vector<Real>>>(
101 : "polygon_origins",
102 : "If 'geometry_type = cylinder', origin(s) about which to curve "
103 : "the polygon corners; if not specified, defaults to (0, 0, 0)");
104 1788 : params.addParam<Real>(
105 : "rotation_angle",
106 1788 : 0,
107 : "If 'geometry_type = cylinder' and when curving corners, the rotation angle (degrees) "
108 : "needed to apply to the original mesh to get a polygon boundary with one side horizontal");
109 :
110 : // TODO: stop-gap solution until the MOOSE reactor module does a better job
111 : // of clearing out lingering sidesets used for stitching, but not for actual BCs/physics
112 1788 : params.addParam<std::vector<BoundaryName>>(
113 : "boundaries_to_rebuild",
114 : "Boundary(s) to retain from the original mesh in the new mesh; if not "
115 : "specified, all original boundaries are kept.");
116 :
117 1788 : params.addParam<bool>(
118 : "retain_original_elem_type",
119 1788 : false,
120 : "Whether to skip the conversion "
121 : "from QUAD9 to QUAD8, or from HEX27 to HEX20, to get into NekRS-compatible element type. "
122 : "This is primarily used to just allow MOOSE's AdvancedExtruderGenerator to extrude Quad9 "
123 : "elements.");
124 :
125 894 : params.addClassDescription(
126 : "Converts MOOSE meshes to element types needed for Nek (Quad8 or Hex20), "
127 : "while optionally preserving curved edges (which were faceted) in the original mesh.");
128 894 : return params;
129 894 : }
130 :
131 447 : NekMeshGenerator::NekMeshGenerator(const InputParameters & params)
132 : : MeshGenerator(params),
133 447 : _input(getMesh("input")),
134 894 : _geometry_type(getParam<MooseEnum>("geometry_type")),
135 894 : _axis(getParam<MooseEnum>("axis")),
136 894 : _curve_corners(getParam<bool>("curve_corners")),
137 894 : _rotation_angle(getParam<Real>("rotation_angle")),
138 894 : _retain_original_elem_type(getParam<bool>("retain_original_elem_type")),
139 1341 : _has_moving_boundary(isParamValid("boundary") || _curve_corners)
140 : {
141 447 : if (_geometry_type == "sphere")
142 : {
143 312 : checkUnusedParam(params,
144 : {"axis",
145 : "curve_corners",
146 : "polygon_sides",
147 : "polygon_size",
148 : "corner_radius",
149 : "polygon_layers",
150 : "polygon_layer_smoothing",
151 : "polygon_boundary",
152 : "polygon_origins",
153 : "rotation_angle"},
154 : "'geometry_type = sphere'");
155 : }
156 : else
157 : {
158 423 : if (_curve_corners)
159 525 : checkRequiredParam(params,
160 : {"polygon_sides", "polygon_size", "polygon_boundary", "corner_radius"},
161 : "'curve_corners' is true");
162 : else
163 3828 : checkUnusedParam(params,
164 : {"polygon_sides",
165 : "polygon_size",
166 : "polygon_boundary",
167 : "corner_radius",
168 : "polygon_layers",
169 : "rotation_angle",
170 : "polygon_layer_smoothing",
171 : "polygon_origins"},
172 : "'curve_corners' is false");
173 : }
174 :
175 894 : if (isParamValid("boundary"))
176 600 : checkRequiredParam(params, "radius", "specifying a 'boundary'");
177 :
178 894 : if (isParamValid("radius"))
179 600 : checkRequiredParam(params, "boundary", "specifying a 'radius'");
180 :
181 1341 : if (isParamValid("radius") || _curve_corners)
182 726 : checkRequiredParam(params, "geometry_type", "specifying a 'radius' or 'curve_corners = true'");
183 :
184 1188 : if (!isParamValid("boundary") && !isParamValid("radius"))
185 1029 : checkUnusedParam(
186 : params, {"axis", "origins", "origins_files", "layers"}, "not setting a 'boundary'");
187 1635 : }
188 :
189 : Point
190 320025 : NekMeshGenerator::projectPoint(const Point & origin, const Point & pt) const
191 : {
192 : Point vec = pt - origin;
193 :
194 320025 : if (_geometry_type == "cylinder")
195 314265 : vec(_axis) = 0.0;
196 :
197 320025 : return vec;
198 : }
199 :
200 : Point
201 115269 : NekMeshGenerator::adjustPointToCircle(const unsigned int & node_id,
202 : Elem * elem,
203 : const Real & radius,
204 : const Point & origin) const
205 : {
206 115269 : auto & node = elem->node_ref(node_id);
207 115269 : const Point pt(node(0), node(1), node(2));
208 :
209 : // project point onto the circle plane and convert to unit vector
210 115269 : Point xy_plane = projectPoint(origin, pt);
211 :
212 : // if the point is exactly on the origin, we don't know in which direction to move
213 : // it so that it lands up on the circle
214 115269 : if (geom_utils::isPointZero(xy_plane))
215 3 : mooseError("Node ID ",
216 : node_id,
217 : " of element ",
218 3 : elem->id(),
219 : " is already on the origin (",
220 : pt(0),
221 : ", ",
222 : pt(1),
223 : ", ",
224 : pt(2),
225 : ").\n"
226 : "This node lacks the nonzero unit vector needed to move it.");
227 :
228 115266 : Point adjustment = xy_plane.unit() * (radius - xy_plane.norm());
229 :
230 : node = pt + adjustment;
231 115266 : return adjustment;
232 : }
233 :
234 : BoundaryID
235 720 : NekMeshGenerator::getBoundaryID(const BoundaryName & name, const MeshBase & mesh) const
236 : {
237 : auto & boundary_info = mesh.get_boundary_info();
238 720 : auto id = MooseMeshUtils::getBoundaryID(name, mesh);
239 :
240 : const auto & all_boundaries = boundary_info.get_boundary_ids();
241 720 : if (all_boundaries.find(id) == all_boundaries.end())
242 15 : mooseError("Boundary '", name, "' was not found in the mesh!");
243 :
244 705 : return id;
245 : }
246 :
247 : void
248 126 : NekMeshGenerator::checkPointLength(const std::vector<std::vector<Real>> & points,
249 : std::string name) const
250 : {
251 306 : for (const auto & o : points)
252 : {
253 186 : if (o.size() == 0)
254 6 : mooseError("Zero-length entry in '" + name +
255 : "' detected! Please be sure that each "
256 3 : "entry in '" +
257 : name + "' has a length\ndivisible by 3 to represent (x, y, z) coordinates.");
258 :
259 183 : if (o.size() % 3 != 0)
260 6 : mooseError("When using multiple origins for one boundary, each entry in '" + name +
261 : "' "
262 : "must have a length\ndivisible by 3 to represent (x, y, z) coordinates!");
263 : }
264 120 : }
265 :
266 : void
267 90513 : NekMeshGenerator::adjustMidPointNode(const unsigned int & node_id, Elem * elem) const
268 : {
269 90513 : std::pair<unsigned int, unsigned int> p = pairedNodesAboutMidPoint(node_id);
270 :
271 90513 : auto & adjust = elem->node_ref(node_id);
272 90513 : const auto & primary = elem->node_ref(p.first);
273 90513 : const auto & secondary = elem->node_ref(p.second);
274 :
275 90513 : Point pt(0.5 * (primary(0) + secondary(0)),
276 90513 : 0.5 * (primary(1) + secondary(1)),
277 90513 : 0.5 * (primary(2) + secondary(2)));
278 :
279 : adjust = pt;
280 90513 : }
281 :
282 : bool
283 45684 : NekMeshGenerator::isNearCorner(const Point & pt) const
284 : {
285 : // point is a moving point by the corner if the distance from the corner is less
286 : // that the distance from the circle tangent to the corner
287 45684 : Real distance_to_closest_corner = geom_utils::minDistanceToPoints(pt, _polygon_corners, _axis);
288 :
289 45684 : return distance_to_closest_corner < _max_corner_distance;
290 : }
291 :
292 : unsigned int
293 146232 : NekMeshGenerator::getNodeIndex(const Elem * elem, const Point & pt) const
294 : {
295 1388232 : for (unsigned int i = 0; i < _n_start_nodes; ++i)
296 1388232 : if (pt.absolute_fuzzy_equals(elem->point(i)))
297 146232 : return i;
298 :
299 0 : mooseError("Failed to find any node on element ", elem->id(), " that matches ", pt);
300 : }
301 :
302 : Point
303 17664 : NekMeshGenerator::getClosestOrigin(const unsigned int & index, const Point & pt) const
304 : {
305 17664 : const auto & candidates = _origin[index];
306 : double distance = std::numeric_limits<double>::max();
307 : int origin_index;
308 :
309 17664 : int n_origins = candidates.size() / 3;
310 222420 : for (int i = 0; i < n_origins; ++i)
311 : {
312 204756 : Point origin(candidates[3 * i], candidates[3 * i + 1], candidates[3 * i + 2]);
313 :
314 : // get the distance to this origin
315 204756 : Point d = projectPoint(origin, pt);
316 204756 : Real current_distance = d.norm();
317 :
318 204756 : if (current_distance < distance)
319 : {
320 : distance = current_distance;
321 : origin_index = i;
322 : }
323 : }
324 :
325 17664 : Point closest(candidates[3 * origin_index],
326 17664 : candidates[3 * origin_index + 1],
327 17664 : candidates[3 * origin_index + 2]);
328 17664 : return closest;
329 : }
330 :
331 : void
332 17670 : NekMeshGenerator::moveElem(Elem * elem,
333 : const unsigned int & boundary_index,
334 : const unsigned int & primary_face,
335 : const std::vector<Real> & polygon_layer_smoothing)
336 : {
337 17670 : bool is_corner_boundary = boundary_index >= _n_noncorner_boundaries;
338 :
339 17670 : const auto bl_elems = getBoundaryLayerElems(elem, _layers[boundary_index], primary_face);
340 :
341 : // use the element centroid for finding the closest origin
342 17664 : const Point centroid = elem->vertex_average();
343 17664 : Point pt = getClosestOrigin(boundary_index, centroid);
344 :
345 163638 : for (auto & face_node : _side_nodes_map[primary_face])
346 : {
347 145977 : const Node & n = elem->node_ref(face_node);
348 145977 : const Point corner_point(n(0), n(1), n(2));
349 145977 : if (is_corner_boundary && !isNearCorner(corner_point))
350 30708 : continue;
351 :
352 : // move the points on the primary face
353 115269 : Point adjustment = adjustPointToCircle(face_node, elem, _radius[boundary_index], pt);
354 :
355 : // move boundary layers of paired nodes, if present
356 : Elem * bl_elem = elem;
357 115266 : unsigned int start_node = face_node;
358 115266 : unsigned int start_face = primary_face;
359 115266 : unsigned int pair_node = pairedFaceNode(start_node, start_face);
360 :
361 261498 : for (unsigned int l = 0; l < _layers[boundary_index]; ++l)
362 : {
363 : auto & paired_node = bl_elem->node_ref(pair_node);
364 146232 : Real multiplier = is_corner_boundary ? polygon_layer_smoothing[l] : 1.0;
365 :
366 : paired_node += adjustment * multiplier;
367 :
368 : // if this is a corner node, we also need to adjust the mid-point node
369 146232 : if (isCornerNode(start_node))
370 68544 : adjustMidPointNode(midPointNodeIndex(start_face, start_node), bl_elem);
371 :
372 : // increment to the next boundary layer element
373 146232 : auto next_elem = bl_elems[l];
374 : bl_elem = const_cast<Elem *>(next_elem);
375 146232 : Point pt(paired_node(0), paired_node(1), paired_node(2));
376 146232 : start_node = getNodeIndex(bl_elem, pt);
377 146232 : pair_node = pairedFaceNode(start_node, start_face);
378 : }
379 :
380 : // even if there aren't boundary layers, we need to adjust the mid-point side node
381 : // of the first moved element
382 115266 : if (_layers[boundary_index] == 0)
383 48954 : if (isCornerNode(face_node))
384 21969 : adjustMidPointNode(midPointNodeIndex(primary_face, face_node), elem);
385 : }
386 17661 : }
387 :
388 : std::vector<Elem *>
389 17670 : NekMeshGenerator::getBoundaryLayerElems(Elem * elem,
390 : const unsigned int & n_layers,
391 : const unsigned int & primary_face) const
392 : {
393 : std::vector<Elem *> nested_elems;
394 :
395 17670 : Elem * bl_elem = elem;
396 17670 : unsigned int start_face = primary_face;
397 17670 : unsigned int pair_face = _across_face[start_face];
398 :
399 38133 : for (unsigned int l = 0; l < n_layers; ++l)
400 : {
401 : // increment to the next boundary layer element
402 20469 : auto next_elem = getNextLayerElem(*bl_elem, pair_face, start_face);
403 20463 : bl_elem = const_cast<Elem *>(next_elem);
404 20463 : pair_face = _across_face[start_face];
405 20463 : nested_elems.push_back(bl_elem);
406 : }
407 :
408 17664 : return nested_elems;
409 : }
410 :
411 : unsigned int
412 261498 : NekMeshGenerator::pairedFaceNode(const unsigned int & node_id, const unsigned int & face_id) const
413 : {
414 : unsigned int pair;
415 :
416 1272891 : for (const auto & p : _across_pair[face_id])
417 : {
418 1272891 : if (p.first == node_id)
419 : {
420 261498 : pair = p.second;
421 261498 : break;
422 : }
423 : }
424 :
425 261498 : return pair;
426 : }
427 :
428 : unsigned int
429 90513 : NekMeshGenerator::midPointNodeIndex(const unsigned int & face_id,
430 : const unsigned int & face_node) const
431 : {
432 90513 : const auto & primary_nodes = _corner_nodes[face_id];
433 90513 : auto it = std::find(primary_nodes.begin(), primary_nodes.end(), face_node);
434 90513 : return _side_ids[face_id][it - primary_nodes.begin()];
435 : }
436 :
437 : const Elem *
438 20469 : NekMeshGenerator::getNextLayerElem(const Elem & elem,
439 : const unsigned int & touching_face,
440 : unsigned int & next_touching_face) const
441 :
442 : {
443 : std::set<const Elem *> neighbor_set;
444 :
445 : // for the input element, get the node index on the mid-face, since for Hex27 that should
446 : // only uniquely touch one other element
447 20469 : auto face_node = getFaceNode(touching_face);
448 20469 : auto face_pt = elem.point(face_node);
449 :
450 20469 : elem.find_point_neighbors(face_pt, neighbor_set);
451 :
452 20469 : if (neighbor_set.size() != 2)
453 6 : mooseError(
454 : "Boundary layer sweeping requires finding exactly one neighbor element\n"
455 : "through the layer face! Found ",
456 6 : neighbor_set.size() - 1,
457 : " neighbors for element ",
458 6 : elem.id(),
459 : ", face ",
460 : touching_face,
461 : ".\n\nThis can happen if you have specified more 'layers' than are actually in your mesh.");
462 :
463 : // When in 3-D, this mesh generator currently assumes we're working on an extruded mesh,
464 : // so that the face pattern follows as we sweep through the boundary layers
465 20463 : next_touching_face = _across_face[touching_face];
466 :
467 : std::set<const Elem *>::iterator it = neighbor_set.begin();
468 38739 : for (std::size_t i = 0; i < neighbor_set.size(); ++i, it++)
469 : {
470 : // we restrict the size to 2, so just return the element that is NOT the input element
471 38739 : if ((*it)->id() != elem.id())
472 20463 : return *it;
473 : }
474 :
475 0 : mooseError("Failed to find a neighbor element across the boundary layer! Please check that\n"
476 : "the 'layers' are set to reasonable values.");
477 : }
478 :
479 : void
480 384 : NekMeshGenerator::moveNodes(std::unique_ptr<MeshBase> & mesh,
481 : std::vector<Real> & polygon_layer_smoothing)
482 : {
483 : auto & boundary_info = mesh->get_boundary_info();
484 :
485 : // move the nodes on the surface of interest
486 166620 : for (auto & elem : mesh->element_ptr_range())
487 : {
488 : bool at_least_one_face_on_boundary = false;
489 :
490 561828 : for (unsigned short int s = 0; s < _n_sides; ++s)
491 : {
492 : // get the boundary IDs that this element face lie on
493 : std::vector<boundary_id_type> b;
494 478902 : boundary_info.boundary_ids(elem, s, b);
495 :
496 : // if there is no moving boundary, or no faces of this element are on any boundaries, we can
497 : // leave
498 478902 : if (!_has_moving_boundary || b.size() == 0)
499 343080 : continue;
500 :
501 : // is this face on any of the specified boundaries?
502 : std::vector<int> indices;
503 :
504 282630 : for (const auto & b_index : b)
505 : {
506 146808 : auto it = std::find(_moving_boundary.begin(), _moving_boundary.end(), b_index);
507 :
508 146808 : if (it != _moving_boundary.end())
509 17685 : indices.push_back(it - _moving_boundary.begin());
510 : }
511 :
512 : // if none of this element's faces are on the boundaries of interest, we can leave
513 135822 : if (indices.size() == 0)
514 : continue;
515 :
516 17682 : if (indices.size() > 1)
517 3 : mooseError("This mesh generator does not support elements with the same face "
518 : "existing on multiple moving side sets!");
519 :
520 17679 : unsigned int index = indices[0];
521 :
522 17679 : if (at_least_one_face_on_boundary)
523 9 : mooseError("This mesh generator cannot be applied to elements that have more than "
524 : "one face on the circular sideset!");
525 :
526 : at_least_one_face_on_boundary = true;
527 17670 : moveElem(elem, index, s, polygon_layer_smoothing);
528 : }
529 363 : }
530 363 : }
531 :
532 : unsigned int
533 20469 : NekMeshGenerator::getFaceNode(const unsigned int & primary_face) const
534 : {
535 20469 : return _side_nodes_map[primary_face][_n_start_nodes_per_side - 1];
536 : }
537 :
538 : void
539 390 : NekMeshGenerator::checkElementType(std::unique_ptr<MeshBase> & mesh)
540 : {
541 : // the type of the first element, should match the type of all other elements
542 390 : _etype = mesh->elem_ptr(0)->type();
543 :
544 169986 : for (const auto & elem : mesh->element_ptr_range())
545 : {
546 84609 : if (_etype != elem->type())
547 0 : mooseError(
548 : "This mesh generator can only be applied to meshes that contain a single element type!");
549 :
550 84609 : switch (elem->type())
551 : {
552 : case QUAD9:
553 : break;
554 : case HEX27:
555 : break;
556 6 : default:
557 6 : mooseError("This mesh generator can only be applied to meshes that contain QUAD9 or HEX27 "
558 : "elements!");
559 : }
560 384 : }
561 384 : }
562 :
563 : void
564 384 : NekMeshGenerator::initializeElemData(std::unique_ptr<MeshBase> & mesh)
565 : {
566 384 : const ElemType etype = mesh->elem_ptr(0)->type();
567 384 : switch (etype)
568 : {
569 144 : case QUAD9:
570 : {
571 144 : _n_start_nodes = Quad9::num_nodes;
572 144 : _n_start_nodes_per_side = Quad9::nodes_per_side;
573 144 : _n_sides = Quad9::num_sides;
574 144 : _n_corner_nodes = Quad4::num_nodes;
575 :
576 144 : if (_retain_original_elem_type)
577 0 : _n_end_nodes = Quad9::num_nodes;
578 : else
579 144 : _n_end_nodes = Quad8::num_nodes;
580 :
581 144 : _side_nodes_map.resize(_n_sides);
582 720 : for (unsigned int i = 0; i < _n_sides; ++i)
583 2304 : for (unsigned int j = 0; j < _n_start_nodes_per_side; ++j)
584 1728 : _side_nodes_map[i].push_back(Quad9::side_nodes_map[i][j]);
585 :
586 144 : _face_nodes_map = _side_nodes_map;
587 :
588 : // for each face, the mid-side nodes to be adjusted
589 144 : _side_ids.push_back({7, 5});
590 144 : _side_ids.push_back({4, 6});
591 144 : _side_ids.push_back({5, 7});
592 144 : _side_ids.push_back({4, 6});
593 :
594 : // corner nodes for each face
595 144 : _corner_nodes.resize(Quad9::num_sides);
596 720 : for (unsigned int i = 0; i < Quad9::num_sides; ++i)
597 1728 : for (unsigned int j = 0; j < Quad4::nodes_per_side; ++j)
598 1152 : _corner_nodes[i].push_back(Quad9::side_nodes_map[i][j]);
599 :
600 144 : _across_pair.resize(Quad9::num_sides);
601 144 : _across_pair[0] = {{0, 3}, {4, 6}, {1, 2}};
602 144 : _across_pair[1] = {{1, 0}, {5, 7}, {2, 3}};
603 144 : _across_pair[2] = {{2, 1}, {6, 4}, {3, 0}};
604 144 : _across_pair[3] = {{3, 2}, {7, 5}, {0, 1}};
605 :
606 144 : _across_face.resize(Quad9::num_sides);
607 144 : _across_face[0] = 2;
608 144 : _across_face[1] = 3;
609 144 : _across_face[2] = 0;
610 144 : _across_face[3] = 1;
611 144 : break;
612 : }
613 240 : case HEX27:
614 : {
615 240 : _n_start_nodes = Hex27::num_nodes;
616 240 : _n_start_nodes_per_side = Hex27::nodes_per_side;
617 240 : _n_sides = Hex27::num_sides;
618 240 : _n_corner_nodes = Hex8::num_nodes;
619 :
620 240 : if (_retain_original_elem_type)
621 0 : _n_end_nodes = Hex27::num_nodes;
622 : else
623 240 : _n_end_nodes = Hex20::num_nodes;
624 :
625 240 : _side_nodes_map.resize(_n_sides);
626 1680 : for (unsigned int i = 0; i < _n_sides; ++i)
627 14400 : for (unsigned int j = 0; j < _n_start_nodes_per_side; ++j)
628 12960 : _side_nodes_map[i].push_back(Hex27::side_nodes_map[i][j]);
629 :
630 240 : _face_nodes_map.resize(Hex27::num_edges);
631 3120 : for (unsigned int i = 0; i < Hex27::num_edges; ++i)
632 11520 : for (unsigned int j = 0; j < Hex27::nodes_per_edge; ++j)
633 8640 : _face_nodes_map[i].push_back(Hex27::edge_nodes_map[i][j]);
634 :
635 : // for each face, the mid-side nodes to be adjusted
636 240 : _side_ids.push_back({12, 15, 14, 13});
637 240 : _side_ids.push_back({11, 9, 17, 19});
638 240 : _side_ids.push_back({8, 18, 10, 16});
639 240 : _side_ids.push_back({9, 11, 19, 17});
640 240 : _side_ids.push_back({10, 8, 16, 18});
641 240 : _side_ids.push_back({12, 13, 14, 15});
642 :
643 : // corner nodes for each face
644 240 : _corner_nodes.resize(Hex27::num_sides);
645 1680 : for (unsigned int i = 0; i < Hex27::num_sides; ++i)
646 7200 : for (unsigned int j = 0; j < Hex8::nodes_per_side; ++j)
647 5760 : _corner_nodes[i].push_back(Hex27::side_nodes_map[i][j]);
648 :
649 240 : _across_pair.resize(Hex27::num_sides);
650 : _across_pair[0] = {
651 240 : {0, 4}, {8, 16}, {1, 5}, {11, 19}, {20, 25}, {9, 17}, {3, 7}, {10, 18}, {2, 6}};
652 : _across_pair[1] = {
653 240 : {0, 3}, {8, 10}, {1, 2}, {12, 15}, {21, 23}, {13, 14}, {4, 7}, {16, 18}, {5, 6}};
654 : _across_pair[2] = {
655 240 : {1, 0}, {9, 11}, {2, 3}, {13, 12}, {22, 24}, {14, 15}, {5, 4}, {17, 19}, {6, 7}};
656 : _across_pair[3] = {
657 240 : {3, 0}, {10, 8}, {2, 1}, {15, 12}, {23, 21}, {14, 13}, {7, 4}, {18, 16}, {6, 5}};
658 : _across_pair[4] = {
659 240 : {0, 1}, {11, 9}, {3, 2}, {12, 13}, {24, 22}, {15, 14}, {4, 5}, {19, 17}, {7, 6}};
660 : _across_pair[5] = {
661 240 : {4, 0}, {16, 8}, {5, 1}, {19, 11}, {25, 20}, {17, 9}, {7, 3}, {18, 10}, {6, 2}};
662 :
663 240 : _across_face.resize(Hex27::num_sides);
664 240 : _across_face[0] = 5;
665 240 : _across_face[1] = 3;
666 240 : _across_face[2] = 4;
667 240 : _across_face[3] = 1;
668 240 : _across_face[4] = 2;
669 240 : _across_face[5] = 0;
670 240 : break;
671 : }
672 0 : default:
673 0 : mooseError("Unhandled element type in initializeElemData()!");
674 : }
675 384 : }
676 :
677 : bool
678 195186 : NekMeshGenerator::isCornerNode(const unsigned int & node) const
679 : {
680 195186 : return node < _n_corner_nodes;
681 : }
682 :
683 : std::pair<unsigned int, unsigned int>
684 90513 : NekMeshGenerator::pairedNodesAboutMidPoint(const unsigned int & node_id) const
685 : {
686 90513 : int index = node_id - _n_corner_nodes;
687 90513 : unsigned int p0 = _face_nodes_map[index][0];
688 90513 : unsigned int p1 = _face_nodes_map[index][1];
689 90513 : return {p0, p1};
690 : }
691 :
692 : std::unique_ptr<MeshBase>
693 447 : NekMeshGenerator::generate()
694 : {
695 447 : std::unique_ptr<MeshBase> mesh = std::move(_input);
696 : auto & boundary_info = mesh->get_boundary_info();
697 : const auto & original_boundaries = boundary_info.get_boundary_ids();
698 :
699 : // TODO: no real reason for this restriction, just didn't need it in the first pass
700 447 : if (!mesh->is_replicated())
701 0 : mooseError("This mesh generator does not yet support distributed mesh implementations!");
702 :
703 : // get the boundary movement information, and check for valid user specifications
704 894 : if (isParamValid("boundary"))
705 : {
706 900 : _radius = getParam<std::vector<Real>>("radius");
707 :
708 672 : for (const auto & r : _radius)
709 375 : if (r <= 0.0)
710 3 : mooseError("All entries in 'radius' must be non-zero and positive!");
711 :
712 594 : const auto & moving_names = getParam<std::vector<BoundaryName>>("boundary");
713 663 : for (const auto & name : moving_names)
714 372 : _moving_boundary.push_back(getBoundaryID(name, *mesh));
715 :
716 291 : if (_moving_boundary.size() != _radius.size())
717 3 : mooseError("'boundary' and 'radius' must be the same length!"
718 : "\n 'boundary' length: ",
719 3 : _moving_boundary.size(),
720 : "\n 'radius' length: ",
721 3 : _radius.size());
722 :
723 798 : if (isParamValid("origins") && isParamValid("origins_files"))
724 3 : mooseError("Cannot specify both 'origins' and 'origins_files'!");
725 :
726 570 : if (isParamValid("origins"))
727 : {
728 324 : _origin = getParam<std::vector<std::vector<Real>>>("origins");
729 :
730 108 : if (_moving_boundary.size() != _origin.size())
731 3 : mooseError("'boundary' and 'origins' must be the same length!"
732 : "\n 'boundary' length: ",
733 3 : _moving_boundary.size(),
734 : "\n 'origins' length: ",
735 3 : _origin.size());
736 :
737 204 : checkPointLength(_origin, "origins");
738 : }
739 354 : else if (isParamValid("origins_files"))
740 : {
741 54 : auto origin_filenames = getParam<std::vector<std::string>>("origins_files");
742 :
743 18 : if (_moving_boundary.size() != origin_filenames.size())
744 3 : mooseError("'boundary' and 'origins_files' must be the same length!"
745 : "\n 'boundary' length: ",
746 3 : _moving_boundary.size(),
747 : "\n 'origins_files' length: ",
748 3 : origin_filenames.size());
749 :
750 15 : _origin.resize(origin_filenames.size());
751 :
752 : int i = 0;
753 27 : for (const auto & f : origin_filenames)
754 : {
755 15 : MooseUtils::DelimitedFileReader file(f, &_communicator);
756 : file.setFormatFlag(MooseUtils::DelimitedFileReader::FormatFlag::ROWS);
757 15 : file.read();
758 :
759 15 : const std::vector<std::vector<double>> & data = file.getData();
760 :
761 99 : for (const auto & d : data)
762 : {
763 87 : if (d.size() != 3)
764 3 : mooseError("All entries in '",
765 : f,
766 : "' must contain exactly 3 entries to represent (x, y, z) coordinates!");
767 :
768 84 : _origin[i].push_back(d[0]);
769 84 : _origin[i].push_back(d[1]);
770 84 : _origin[i].push_back(d[2]);
771 : }
772 :
773 12 : i += 1;
774 12 : }
775 12 : }
776 : else
777 : {
778 : // set to the default value of (0, 0, 0)
779 354 : for (std::size_t i = 0; i < _radius.size(); ++i)
780 390 : _origin.push_back({0.0, 0.0, 0.0});
781 : }
782 :
783 540 : if (isParamValid("layers"))
784 : {
785 180 : _layers = getParam<std::vector<unsigned int>>("layers");
786 :
787 60 : if (_moving_boundary.size() != _layers.size())
788 3 : mooseError("'boundary' and 'layers' must be the same length!"
789 : "\n 'boundary' length: ",
790 3 : _moving_boundary.size(),
791 : "\n 'layers' length: ",
792 3 : _layers.size());
793 : }
794 : else
795 : {
796 : // set to the default values of 0
797 480 : for (std::size_t i = 0; i < _moving_boundary.size(); ++i)
798 270 : _layers.push_back(0.0);
799 : }
800 : }
801 :
802 : // get information related to moving corners
803 : std::vector<Real> polygon_layer_smoothing;
804 414 : _n_noncorner_boundaries = _moving_boundary.size();
805 :
806 414 : if (_curve_corners)
807 : {
808 150 : auto polygon_sides = getParam<unsigned int>("polygon_sides");
809 150 : auto polygon_size = getParam<Real>("polygon_size");
810 150 : auto corner_radius = getParam<Real>("corner_radius");
811 150 : auto polygon_layers = getParam<unsigned int>("polygon_layers");
812 :
813 : std::vector<std::vector<Real>> polygon_origin;
814 150 : if (isParamValid("polygon_origins"))
815 : {
816 42 : polygon_origin = getParam<std::vector<std::vector<Real>>>("polygon_origins");
817 21 : checkPointLength(polygon_origin, "polygon_origins");
818 :
819 : // Cannot specify a non-zero rotation when giving a different polygon origin,
820 : // because we simply rotate the point about either the x, y, or z 'axis'. This
821 : // is not a hard limit, just something we didn't initially put effort to support.
822 21 : if (std::abs(_rotation_angle) > 1e-8)
823 3 : mooseError(
824 : "Cannot specify a non-zero 'rotation_angle' when providing custom 'polygon_origins'!");
825 : }
826 : else
827 108 : polygon_origin.push_back({0.0, 0.0, 0.0});
828 :
829 72 : if (polygon_layers)
830 : {
831 30 : if (isParamValid("polygon_layer_smoothing"))
832 : {
833 30 : polygon_layer_smoothing = getParam<std::vector<Real>>("polygon_layer_smoothing");
834 15 : if (polygon_layers != polygon_layer_smoothing.size())
835 3 : mooseError("The length of 'polygon_layer_smoothing' must be equal to 'polygon_layers'!");
836 :
837 42 : for (auto & p : polygon_layer_smoothing)
838 33 : if (p <= 0.0)
839 3 : mooseError("Each entry in 'polygon_layer_smoothing' must be positive and non-zero!");
840 : }
841 : else
842 : {
843 0 : for (unsigned int i = 0; i < polygon_layers; ++i)
844 0 : polygon_layer_smoothing.push_back(1.0);
845 : }
846 : }
847 : else
848 114 : checkUnusedParam(parameters(), "polygon_layer_smoothing", "not setting 'polygon_layers'");
849 :
850 66 : Real polygon_angle = M_PI - (2.0 * M_PI / polygon_sides);
851 66 : Real max_circle_radius = polygon_size * std::cos(M_PI / polygon_sides);
852 :
853 66 : if (corner_radius > max_circle_radius)
854 3 : mooseError("Specified 'corner_radius' cannot fit within the specified polygon!\n"
855 : "The maximum allowable radius of curvature is: ",
856 : max_circle_radius);
857 :
858 126 : const auto & name = getParam<BoundaryName>("polygon_boundary");
859 63 : auto polygon_boundary = getBoundaryID(name, *mesh);
860 :
861 : // polygon boundary shouldn't already have been specified in the 'boundary'
862 60 : if (std::count(_moving_boundary.begin(), _moving_boundary.end(), polygon_boundary))
863 3 : mooseError("The 'polygon_boundary' cannot also be listed in the 'boundary'!");
864 :
865 57 : Real theta = M_PI / 2.0 - polygon_angle / 2.0;
866 57 : Real l = corner_radius / std::cos(theta);
867 57 : _max_corner_distance = l * std::sin(theta);
868 :
869 : // find origins of the cylinders for the corner fitting
870 : std::vector<Point> corner_origins;
871 150 : for (const auto & o : polygon_origin)
872 : {
873 93 : Point shift(o[0], o[1], o[2]);
874 93 : auto tmp1 = geom_utils::polygonCorners(polygon_sides, polygon_size, _axis);
875 651 : for (const auto & t : tmp1)
876 558 : _polygon_corners.push_back(t + shift);
877 :
878 93 : auto tmp2 = geom_utils::polygonCorners(polygon_sides, polygon_size - l, _axis);
879 651 : for (const auto & t : tmp2)
880 558 : corner_origins.push_back(t + shift);
881 : }
882 : // apply optional rotation
883 57 : Real rotation_angle_radians = _rotation_angle * M_PI / 180.0;
884 : Point axis(0.0, 0.0, 0.0);
885 57 : axis(_axis) = 1.0;
886 615 : for (auto & o : _polygon_corners)
887 558 : o = geom_utils::rotatePointAboutAxis(o, rotation_angle_radians, axis);
888 615 : for (auto & o : corner_origins)
889 558 : o = geom_utils::rotatePointAboutAxis(o, rotation_angle_radians, axis);
890 :
891 : std::vector<Real> flattened_corner_origins;
892 615 : for (const auto & o : corner_origins)
893 2232 : for (unsigned int i = 0; i < 3; ++i)
894 1674 : flattened_corner_origins.push_back(o(i));
895 :
896 : // We can treat the polygon corners simply as extra entries in the
897 : // boundary, origins, and radii vectors
898 57 : _moving_boundary.push_back(polygon_boundary);
899 57 : _radius.push_back(corner_radius);
900 57 : _origin.push_back(flattened_corner_origins);
901 171 : _layers.push_back(getParam<unsigned int>("polygon_layers"));
902 57 : }
903 :
904 : // get information on which boundaries to rebuild, and check for valid user specifications
905 792 : if (isParamValid("boundaries_to_rebuild"))
906 : {
907 234 : const auto & rebuild_names = getParam<std::vector<BoundaryName>>("boundaries_to_rebuild");
908 396 : for (const auto & name : rebuild_names)
909 285 : _boundaries_to_rebuild.insert(getBoundaryID(name, *mesh));
910 : }
911 : else
912 : {
913 : // by default, rebuild all the boundaries
914 2112 : for (const auto & b : original_boundaries)
915 1833 : _boundaries_to_rebuild.insert(b);
916 : }
917 :
918 : // check for valid element type
919 390 : checkElementType(mesh);
920 :
921 384 : initializeElemData(mesh);
922 :
923 : // store all information from the incoming mesh that is needed to rebuild it from scratch
924 : std::vector<dof_id_type> boundary_elem_ids;
925 : std::vector<unsigned int> boundary_face_ids;
926 : std::vector<std::vector<boundary_id_type>> boundary_ids;
927 : std::vector<std::vector<dof_id_type>> node_ids;
928 : std::vector<dof_id_type> elem_ids;
929 : std::vector<subdomain_id_type> elem_block_ids;
930 384 : node_ids.resize(mesh->n_elem());
931 :
932 : // store boundary ID <-> name mapping
933 3357 : for (const auto & b : original_boundaries)
934 8919 : _boundary_id_to_name[b] = boundary_info.get_sideset_name(b);
935 :
936 169974 : for (auto & elem : mesh->element_ptr_range())
937 : {
938 : // store information about the element faces
939 571683 : for (unsigned short int s = 0; s < _n_sides; ++s)
940 : {
941 : std::vector<boundary_id_type> b;
942 487080 : boundary_info.boundary_ids(elem, s, b);
943 :
944 487080 : boundary_elem_ids.push_back(elem->id());
945 487080 : boundary_face_ids.push_back(s);
946 487080 : boundary_ids.push_back(b);
947 : }
948 384 : }
949 :
950 : int i = 0;
951 85371 : for (auto & elem : mesh->element_ptr_range())
952 : {
953 : // store information about the elements
954 84603 : elem_ids.push_back(elem->id());
955 84603 : elem_block_ids.push_back(elem->subdomain_id());
956 :
957 1653435 : for (unsigned int j = 0; j < _n_end_nodes; ++j)
958 1568832 : node_ids[i].push_back(elem->node_ref(j).id());
959 :
960 84603 : i++;
961 384 : }
962 :
963 : // move the nodes on the surface of interest
964 384 : moveNodes(mesh, polygon_layer_smoothing);
965 :
966 : // save and rebuild the nodes
967 : std::vector<Node> original_nodes;
968 : std::vector<dof_id_type> original_node_ids;
969 :
970 : // store information about the nodes
971 795630 : for (const auto & node : mesh->node_ptr_range())
972 : {
973 794904 : original_nodes.push_back(*node);
974 794904 : original_node_ids.push_back(node->id());
975 363 : }
976 :
977 363 : mesh->clear();
978 :
979 : // create the nodes
980 795267 : for (unsigned int i = 0; i < original_nodes.size(); ++i)
981 : {
982 794904 : Point pt(original_nodes[i](0), original_nodes[i](1), original_nodes[i](2));
983 794904 : mesh->add_point(pt, original_node_ids[i]);
984 : }
985 :
986 : // create the elements
987 83247 : for (unsigned int i = 0; i < elem_ids.size(); ++i)
988 : {
989 : Elem * elem;
990 82884 : switch (_etype)
991 : {
992 9324 : case QUAD9:
993 : {
994 9324 : if (_retain_original_elem_type)
995 0 : elem = new Quad9;
996 : else
997 9324 : elem = new Quad8;
998 : break;
999 : }
1000 73560 : case HEX27:
1001 : {
1002 73560 : if (_retain_original_elem_type)
1003 0 : elem = new Hex27;
1004 : else
1005 73560 : elem = new Hex20;
1006 : break;
1007 : }
1008 0 : default:
1009 0 : mooseError("Unhandled element type in generate()!");
1010 : }
1011 :
1012 82884 : elem->set_id(elem_ids[i]);
1013 82884 : elem->subdomain_id() = elem_block_ids[i];
1014 82884 : mesh->add_elem(elem);
1015 :
1016 : const auto & ids = node_ids[i];
1017 1628676 : for (unsigned int n = 0; n < ids.size(); ++n)
1018 : {
1019 1545792 : auto node_ptr = mesh->node_ptr(ids[n]);
1020 1545792 : elem->set_node(n) = node_ptr;
1021 : }
1022 : }
1023 :
1024 : // create the sidesets
1025 479019 : for (unsigned int i = 0; i < boundary_elem_ids.size(); ++i)
1026 : {
1027 478656 : auto elem_id = boundary_elem_ids[i];
1028 478656 : auto boundary = boundary_ids[i];
1029 478656 : const auto & elem = mesh->elem_ptr(elem_id);
1030 :
1031 : // if boundary is not included in the list to rebuild, skip it. Otherwise,
1032 : // rebuild it
1033 629376 : for (const auto & b : boundary)
1034 : {
1035 150720 : if (_boundaries_to_rebuild.find(b) == _boundaries_to_rebuild.end())
1036 59184 : continue;
1037 : else
1038 91536 : boundary_info.add_side(elem, boundary_face_ids[i], b);
1039 : }
1040 : }
1041 :
1042 3153 : for (const auto & b : _boundary_id_to_name)
1043 : {
1044 : // if boundary is not included in the list to rebuild, skip it
1045 2790 : if (_boundaries_to_rebuild.find(b.first) == _boundaries_to_rebuild.end())
1046 867 : continue;
1047 :
1048 1923 : boundary_info.sideset_name(b.first) = b.second;
1049 : }
1050 :
1051 363 : mesh->prepare_for_use();
1052 :
1053 726 : return dynamic_pointer_cast<MeshBase>(mesh);
1054 1089 : }
|
