11 #include "libmesh/face_quad4.h" 13 #include "libmesh/mesh_modification.h" 14 #include "libmesh/serial_mesh.h" 15 #include "libmesh/boundary_info.h" 16 #include "libmesh/utility.h" 27 "full top_right top_left bottom_left bottom_right right_half left_half top_half bottom_half",
30 "num_sectors % 2 = 0, num_sectors > 0" 31 "Number of azimuthal sectors in each quadrant" 32 "'num_sectors' must be an even number.");
33 params.
addRequiredParam<std::vector<Real>>(
"radii",
"Radii of major concentric circles");
35 "rings",
"Number of rings in each circle or in the moderator");
37 "Length of inner square / radius of the innermost circle");
40 "It determines if meshes for a outer square are added to concentric circle meshes.");
45 "The moderator can be added to complete meshes for one unit cell of fuel assembly." 46 "Elements are quad meshes.");
47 params.
addParam<
MooseEnum>(
"portion", portion,
"Control of which part of mesh is created");
49 "preserve_volumes",
"Volume of concentric circles can be preserved using this function.");
50 params.
addClassDescription(
"This ConcentricCircleMesh source code is to generate concentric " 57 _num_sectors(getParam<unsigned
int>(
"num_sectors")),
58 _radii(getParam<
std::vector<
Real>>(
"radii")),
59 _rings(getParam<
std::vector<unsigned
int>>(
"rings")),
60 _inner_mesh_fraction(getParam<
Real>(
"inner_mesh_fraction")),
61 _has_outer_square(getParam<bool>(
"has_outer_square")),
62 _pitch(getParam<
Real>(
"pitch")),
63 _preserve_volumes(getParam<bool>(
"preserve_volumes")),
68 mooseError(
"ConcentricCircleMesh: num_sectors must be an even number.");
71 for (
unsigned i = 0; i <
_radii.size() - 1; ++i)
73 mooseError(
"Radii must be provided in order by starting with the smallest radius and " 74 "providing the following gradual radii.");
78 mooseError(
"The aspect ratio can not be larger than cos(PI/4).");
84 mooseError(
"The size of 'rings' must be equal to the size of 'radii' plus 1.");
89 mooseError(
"The size of 'rings' must be equal to the size of 'radii'.");
93 for (
unsigned i = 0; i <
_radii.size(); ++i)
95 mooseError(
"The pitch / 2 must be larger than any radii.");
102 ReplicatedMesh &
mesh = cast_ref<ReplicatedMesh &>(
getMesh());
104 mesh.set_mesh_dimension(2);
105 mesh.set_spatial_dimension(2);
106 BoundaryInfo & boundary_info =
mesh.get_boundary_info();
110 std::vector<Real> total_concentric_circles;
126 total_concentric_circles.push_back(
_radii[j - 1] +
138 Real original_radius = 0.0;
139 for (
unsigned i = 0; i < total_concentric_circles.size(); ++i)
146 original_radius = total_concentric_circles[i];
147 total_concentric_circles[i] = modified_radius;
156 original_radius = total_concentric_circles[i];
157 total_concentric_circles[i] = modified_radius;
163 unsigned num_total_nodes = 0;
172 std::vector<Node *> nodes(num_total_nodes);
173 unsigned node_id = 0;
182 nodes[node_id] =
mesh.add_point(Point(x, y, 0.0), node_id);
188 Real current_radius = 0.0;
190 for (
unsigned layers = 0; layers < total_concentric_circles.size(); ++layers)
192 current_radius = total_concentric_circles[layers];
193 for (
unsigned num_outer_nodes = 0; num_outer_nodes <=
_num_sectors; ++num_outer_nodes)
195 const Real x = current_radius *
std::cos(num_outer_nodes * d_angle);
196 const Real y = current_radius *
std::sin(num_outer_nodes * d_angle);
197 nodes[node_id] =
mesh.add_point(Point(x, y, 0.0), node_id);
205 Real current_radius_moderator = 0.0;
206 for (
unsigned i = 1; i <=
_rings.back(); ++i)
208 current_radius_moderator =
210 total_concentric_circles.push_back(current_radius_moderator);
211 for (
unsigned num_outer_nodes = 0; num_outer_nodes <=
_num_sectors; ++num_outer_nodes)
213 const Real x = current_radius_moderator *
std::cos(num_outer_nodes * d_angle);
214 const Real y = current_radius_moderator *
std::sin(num_outer_nodes * d_angle);
215 nodes[node_id] =
mesh.add_point(Point(x, y, 0.0), node_id);
223 const Real y =
_pitch / 2 * std::tan(j * d_angle);
224 nodes[node_id] =
mesh.add_point(Point(x, y, 0.0), node_id);
233 nodes[node_id] =
mesh.add_point(Point(x, y, 0.0), node_id);
246 int additional_term = 0;
247 int counter = standard;
250 counter = counter - 2;
251 additional_term = additional_term + counter;
253 limit = standard + additional_term;
255 else if (standard == 4)
259 std::vector<unsigned int> subdomainIDs;
260 for (
unsigned int i = 0; i <
_rings.size(); ++i)
261 for (
unsigned int j = 0; j <
_rings[i]; ++j)
262 subdomainIDs.push_back(i + 1);
265 subdomainIDs.push_back(subdomainIDs.back());
267 while (index <= limit)
269 Elem *
elem =
mesh.add_elem(
new Quad4);
270 elem->set_node(0) = nodes[index];
273 elem->set_node(3) = nodes[index + 1];
274 elem->subdomain_id() = subdomainIDs[0];
276 if (index < standard / 2)
277 boundary_info.add_side(
elem, 3, 1);
278 if (index % (standard / 2 + 1) == 0)
279 boundary_info.add_side(
elem, 0, 2);
282 if ((index - standard / 2) % (standard / 2 + 1) == 0)
290 while (index < limit)
292 Elem *
elem =
mesh.add_elem(
new Quad4);
293 elem->set_node(0) = nodes[index];
296 elem->set_node(3) = nodes[index + 1];
297 elem->subdomain_id() = subdomainIDs[0];
299 if (index == (standard / 2 + 1) * (standard / 2))
300 boundary_info.add_side(
elem, 0, 2);
307 while (index != standard / 2)
309 Elem *
elem =
mesh.add_elem(
new Quad4);
310 elem->set_node(0) = nodes[index];
315 elem->subdomain_id() = subdomainIDs[0];
317 if (index == standard + 1)
318 boundary_info.add_side(
elem, 2, 1);
326 limit =
static_cast<int>(num_total_nodes) - standard - 2;
330 while (index < limit)
333 Elem *
elem =
mesh.add_elem(
new Quad4);
334 elem->set_node(0) = nodes[index];
337 elem->set_node(3) = nodes[index + 1];
339 for (
int i = 0; i < static_cast<int>(subdomainIDs.size()) - 1; ++i)
340 if (index < limit - (standard + 1) * i && index >= limit - (standard + 1) * (i + 1))
341 elem->subdomain_id() = subdomainIDs[subdomainIDs.size() - 1 - i];
346 if ((index - initial) % (standard + 1) == 0)
347 boundary_info.add_side(
elem, 0, 2);
348 if ((index -
final) % (standard + 1) == 0)
349 boundary_info.add_side(
elem, 2, 1);
350 if (index > limit - (standard + 1))
354 if (index < limit - standard + standard / 2)
355 boundary_info.add_side(
elem, 1, 3);
357 boundary_info.add_side(
elem, 1, 4);
361 boundary_info.add_side(
elem, 1, 3);
365 if (index == (num_nodes_boundary + counter * (standard + 1)) - 1)
373 boundary_info.sideset_name(1) =
"left";
374 boundary_info.sideset_name(2) =
"bottom";
377 boundary_info.sideset_name(3) =
"outer";
380 boundary_info.sideset_name(3) =
"right";
381 boundary_info.sideset_name(4) =
"top";
386 MeshTools::Modification::rotate(
mesh, 90, 0, 0);
387 boundary_info.sideset_name(1) =
"bottom";
388 boundary_info.sideset_name(2) =
"right";
391 boundary_info.sideset_name(3) =
"outer";
394 boundary_info.sideset_name(3) =
"top";
395 boundary_info.sideset_name(4) =
"left";
400 MeshTools::Modification::rotate(
mesh, 180, 0, 0);
401 boundary_info.sideset_name(1) =
"right";
402 boundary_info.sideset_name(2) =
"top";
405 boundary_info.sideset_name(3) =
"outer";
408 boundary_info.sideset_name(3) =
"left";
409 boundary_info.sideset_name(4) =
"bottom";
412 else if (
_portion ==
"bottom_right")
414 MeshTools::Modification::rotate(
mesh, 270, 0, 0);
415 boundary_info.sideset_name(1) =
"top";
416 boundary_info.sideset_name(2) =
"left";
419 boundary_info.sideset_name(3) =
"outer";
422 boundary_info.sideset_name(3) =
"bottom";
423 boundary_info.sideset_name(4) =
"right";
429 ReplicatedMesh other_mesh(
mesh);
431 MeshTools::Modification::rotate(other_mesh, 90, 0, 0);
434 MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
435 MeshTools::Modification::change_boundary_id(other_mesh, 2, 6);
436 MeshTools::Modification::change_boundary_id(other_mesh, 3, 7);
437 MeshTools::Modification::change_boundary_id(other_mesh, 4, 1);
438 MeshTools::Modification::change_boundary_id(other_mesh, 5, 2);
439 MeshTools::Modification::change_boundary_id(other_mesh, 6, 3);
440 MeshTools::Modification::change_boundary_id(other_mesh, 7, 4);
441 mesh.prepare_for_use();
442 other_mesh.prepare_for_use();
443 mesh.stitch_meshes(other_mesh, 1, 3, TOLERANCE,
true);
444 mesh.get_boundary_info().sideset_name(1) =
"left";
445 mesh.get_boundary_info().sideset_name(2) =
"bottom";
446 mesh.get_boundary_info().sideset_name(3) =
"right";
447 mesh.get_boundary_info().sideset_name(4) =
"top";
451 MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
452 MeshTools::Modification::change_boundary_id(other_mesh, 2, 1);
453 MeshTools::Modification::change_boundary_id(other_mesh, 5, 2);
454 mesh.prepare_for_use();
455 other_mesh.prepare_for_use();
456 mesh.stitch_meshes(other_mesh, 1, 1, TOLERANCE,
true);
458 MeshTools::Modification::change_boundary_id(
mesh, 2, 1);
459 MeshTools::Modification::change_boundary_id(
mesh, 3, 2);
460 mesh.get_boundary_info().sideset_name(1) =
"bottom";
461 mesh.get_boundary_info().sideset_name(2) =
"outer";
468 ReplicatedMesh other_mesh(
mesh);
470 MeshTools::Modification::rotate(other_mesh, 270, 0, 0);
473 MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
474 MeshTools::Modification::change_boundary_id(other_mesh, 2, 6);
475 MeshTools::Modification::change_boundary_id(other_mesh, 3, 7);
476 MeshTools::Modification::change_boundary_id(other_mesh, 4, 3);
477 MeshTools::Modification::change_boundary_id(other_mesh, 5, 4);
478 MeshTools::Modification::change_boundary_id(other_mesh, 6, 1);
479 MeshTools::Modification::change_boundary_id(other_mesh, 7, 2);
480 mesh.prepare_for_use();
481 other_mesh.prepare_for_use();
482 mesh.stitch_meshes(other_mesh, 2, 4, TOLERANCE,
true);
483 mesh.get_boundary_info().sideset_name(1) =
"left";
484 mesh.get_boundary_info().sideset_name(2) =
"bottom";
485 mesh.get_boundary_info().sideset_name(3) =
"right";
486 mesh.get_boundary_info().sideset_name(4) =
"top";
490 MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
491 MeshTools::Modification::change_boundary_id(other_mesh, 2, 1);
492 MeshTools::Modification::change_boundary_id(other_mesh, 5, 2);
493 mesh.prepare_for_use();
494 other_mesh.prepare_for_use();
495 mesh.stitch_meshes(other_mesh, 2, 2, TOLERANCE,
true);
497 MeshTools::Modification::change_boundary_id(
mesh, 3, 2);
498 mesh.get_boundary_info().sideset_name(1) =
"left";
499 mesh.get_boundary_info().sideset_name(2) =
"outer";
505 ReplicatedMesh other_mesh(
mesh);
508 MeshTools::Modification::rotate(other_mesh, 90, 0, 0);
509 MeshTools::Modification::rotate(
mesh, 180, 0, 0);
513 MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
514 MeshTools::Modification::change_boundary_id(other_mesh, 2, 6);
515 MeshTools::Modification::change_boundary_id(other_mesh, 3, 7);
516 MeshTools::Modification::change_boundary_id(other_mesh, 4, 1);
517 MeshTools::Modification::change_boundary_id(other_mesh, 5, 2);
518 MeshTools::Modification::change_boundary_id(other_mesh, 6, 3);
519 MeshTools::Modification::change_boundary_id(other_mesh, 7, 4);
521 MeshTools::Modification::change_boundary_id(
mesh, 1, 5);
522 MeshTools::Modification::change_boundary_id(
mesh, 2, 6);
523 MeshTools::Modification::change_boundary_id(
mesh, 3, 7);
524 MeshTools::Modification::change_boundary_id(
mesh, 4, 2);
525 MeshTools::Modification::change_boundary_id(
mesh, 5, 3);
526 MeshTools::Modification::change_boundary_id(
mesh, 6, 4);
527 MeshTools::Modification::change_boundary_id(
mesh, 7, 1);
528 mesh.prepare_for_use();
529 other_mesh.prepare_for_use();
530 mesh.stitch_meshes(other_mesh, 4, 2, TOLERANCE,
true);
531 mesh.get_boundary_info().sideset_name(1) =
"left";
532 mesh.get_boundary_info().sideset_name(2) =
"bottom";
533 mesh.get_boundary_info().sideset_name(3) =
"right";
534 mesh.get_boundary_info().sideset_name(4) =
"top";
538 MeshTools::Modification::change_boundary_id(
mesh, 1, 5);
539 MeshTools::Modification::change_boundary_id(
mesh, 2, 1);
540 MeshTools::Modification::change_boundary_id(
mesh, 5, 2);
541 mesh.prepare_for_use();
542 other_mesh.prepare_for_use();
543 mesh.stitch_meshes(other_mesh, 1, 1, TOLERANCE,
true);
545 MeshTools::Modification::change_boundary_id(
mesh, 2, 1);
546 MeshTools::Modification::change_boundary_id(
mesh, 3, 2);
547 mesh.get_boundary_info().sideset_name(1) =
"right";
548 mesh.get_boundary_info().sideset_name(2) =
"outer";
554 ReplicatedMesh other_mesh(
mesh);
556 MeshTools::Modification::rotate(other_mesh, 180, 0, 0);
557 MeshTools::Modification::rotate(
mesh, 270, 0, 0);
561 MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
562 MeshTools::Modification::change_boundary_id(other_mesh, 2, 6);
563 MeshTools::Modification::change_boundary_id(other_mesh, 3, 7);
564 MeshTools::Modification::change_boundary_id(other_mesh, 4, 2);
565 MeshTools::Modification::change_boundary_id(other_mesh, 5, 3);
566 MeshTools::Modification::change_boundary_id(other_mesh, 6, 4);
567 MeshTools::Modification::change_boundary_id(other_mesh, 7, 1);
569 MeshTools::Modification::change_boundary_id(
mesh, 1, 5);
570 MeshTools::Modification::change_boundary_id(
mesh, 2, 6);
571 MeshTools::Modification::change_boundary_id(
mesh, 3, 7);
572 MeshTools::Modification::change_boundary_id(
mesh, 4, 3);
573 MeshTools::Modification::change_boundary_id(
mesh, 5, 4);
574 MeshTools::Modification::change_boundary_id(
mesh, 6, 1);
575 MeshTools::Modification::change_boundary_id(
mesh, 7, 2);
576 mesh.prepare_for_use();
577 other_mesh.prepare_for_use();
578 mesh.stitch_meshes(other_mesh, 1, 3, TOLERANCE,
true);
579 mesh.get_boundary_info().sideset_name(1) =
"left";
580 mesh.get_boundary_info().sideset_name(2) =
"bottom";
581 mesh.get_boundary_info().sideset_name(3) =
"right";
582 mesh.get_boundary_info().sideset_name(4) =
"top";
586 MeshTools::Modification::change_boundary_id(
mesh, 1, 5);
587 MeshTools::Modification::change_boundary_id(
mesh, 2, 1);
588 MeshTools::Modification::change_boundary_id(
mesh, 5, 2);
589 mesh.prepare_for_use();
590 other_mesh.prepare_for_use();
591 mesh.stitch_meshes(other_mesh, 1, 1, TOLERANCE,
true);
593 MeshTools::Modification::change_boundary_id(
mesh, 2, 1);
594 MeshTools::Modification::change_boundary_id(
mesh, 3, 2);
595 mesh.get_boundary_info().sideset_name(1) =
"top";
596 mesh.get_boundary_info().sideset_name(2) =
"outer";
602 ReplicatedMesh portion_two(
mesh);
605 MeshTools::Modification::rotate(portion_two, 90, 0, 0);
610 MeshTools::Modification::change_boundary_id(portion_two, 1, 5);
611 MeshTools::Modification::change_boundary_id(portion_two, 2, 6);
612 MeshTools::Modification::change_boundary_id(portion_two, 3, 7);
613 MeshTools::Modification::change_boundary_id(portion_two, 4, 1);
614 MeshTools::Modification::change_boundary_id(portion_two, 5, 2);
615 MeshTools::Modification::change_boundary_id(portion_two, 6, 3);
616 MeshTools::Modification::change_boundary_id(portion_two, 7, 4);
617 mesh.prepare_for_use();
618 portion_two.prepare_for_use();
620 mesh.stitch_meshes(portion_two, 1, 3, TOLERANCE,
true);
623 ReplicatedMesh portion_bottom(
mesh);
624 MeshTools::Modification::rotate(portion_bottom, 180, 0, 0);
625 MeshTools::Modification::change_boundary_id(portion_bottom, 1, 5);
626 MeshTools::Modification::change_boundary_id(portion_bottom, 2, 6);
627 MeshTools::Modification::change_boundary_id(portion_bottom, 3, 7);
628 MeshTools::Modification::change_boundary_id(portion_bottom, 4, 2);
629 MeshTools::Modification::change_boundary_id(portion_bottom, 5, 3);
630 MeshTools::Modification::change_boundary_id(portion_bottom, 6, 4);
631 MeshTools::Modification::change_boundary_id(portion_bottom, 7, 1);
632 mesh.prepare_for_use();
633 portion_bottom.prepare_for_use();
635 mesh.stitch_meshes(portion_bottom, 2, 4, TOLERANCE,
true);
637 mesh.get_boundary_info().sideset_name(1) =
"left";
638 mesh.get_boundary_info().sideset_name(2) =
"bottom";
639 mesh.get_boundary_info().sideset_name(3) =
"right";
640 mesh.get_boundary_info().sideset_name(4) =
"top";
641 portion_bottom.clear();
645 MeshTools::Modification::change_boundary_id(portion_two, 1, 5);
646 MeshTools::Modification::change_boundary_id(portion_two, 2, 1);
647 MeshTools::Modification::change_boundary_id(portion_two, 5, 2);
649 mesh.prepare_for_use();
650 portion_two.prepare_for_use();
651 mesh.stitch_meshes(portion_two, 1, 1, TOLERANCE,
true);
653 ReplicatedMesh portion_bottom(
mesh);
654 MeshTools::Modification::rotate(portion_bottom, 180, 0, 0);
656 mesh.prepare_for_use();
657 portion_bottom.prepare_for_use();
658 mesh.stitch_meshes(portion_bottom, 2, 2, TOLERANCE,
true);
659 MeshTools::Modification::change_boundary_id(
mesh, 3, 1);
660 mesh.get_boundary_info().sideset_name(1) =
"outer";
661 portion_bottom.clear();
667 mesh.prepare_for_use();
static InputParameters validParams()
Typical "Moose-style" constructor and copy constructor.
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template * sin(_arg) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(tan
Real _inner_mesh_fraction
Size of inner square in relation to radius of the innermost concentric circle.
virtual void buildMesh() override
Must be overridden by child classes.
static InputParameters validParams()
unsigned int _num_sectors
Number of sectors in one quadrant.
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
std::vector< Real > _radii
Radii of concentric circles.
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template cos(_arg) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(cos
MeshBase & getMesh()
Accessor for the underlying libMesh Mesh object.
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
ConcentricCircleMesh(const InputParameters ¶meters)
bool _has_outer_square
Adding the moderator is optional.
std::vector< unsigned int > _rings
Number of rings in each circle or in the moderator.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
registerMooseObject("MooseApp", ConcentricCircleMesh)
void mooseError(Args &&... args) const
Emits an error prefixed with object name and type.
Mesh generated from parameters.
virtual Elem * elem(const dof_id_type i)
Various accessors (pointers/references) for Elem "i".
bool _preserve_volumes
Volume preserving function is optional.
void ErrorVector unsigned int
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template pow< 2 >(tan(_arg))+1.0) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(sqrt
MooseEnum _portion
Control of which portion of mesh will be developed.