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 "ConcentricCircleMesh.h"
11 : #include "libmesh/face_quad4.h"
12 : #include "MooseMesh.h"
13 : #include "MooseApp.h"
14 : #include "libmesh/mesh_modification.h"
15 : #include "libmesh/serial_mesh.h"
16 : #include "libmesh/boundary_info.h"
17 : #include "libmesh/utility.h"
18 : // C++ includes
19 : #include <cmath> // provides round, not std::round (see http://www.cplusplus.com/reference/cmath/round/)
20 :
21 : registerMooseObject("MooseApp", ConcentricCircleMesh);
22 :
23 : InputParameters
24 14285 : ConcentricCircleMesh::validParams()
25 : {
26 14285 : InputParameters params = MooseMesh::validParams();
27 : MooseEnum portion(
28 : "full top_right top_left bottom_left bottom_right right_half left_half top_half bottom_half",
29 14285 : "full");
30 14285 : params.addRequiredParam<unsigned int>("num_sectors",
31 : "num_sectors % 2 = 0, num_sectors > 0"
32 : "Number of azimuthal sectors in each quadrant"
33 : "'num_sectors' must be an even number.");
34 14285 : params.addRequiredParam<std::vector<Real>>("radii", "Radii of major concentric circles");
35 14285 : params.addRequiredParam<std::vector<unsigned int>>(
36 : "rings", "Number of rings in each circle or in the moderator");
37 14285 : params.addRequiredParam<Real>("inner_mesh_fraction",
38 : "Length of inner square / radius of the innermost circle");
39 14285 : params.addRequiredParam<bool>(
40 : "has_outer_square",
41 : "It determines if meshes for a outer square are added to concentric circle meshes.");
42 42855 : params.addRangeCheckedParam<Real>(
43 : "pitch",
44 28570 : 0.0,
45 : "pitch>=0.0",
46 : "The moderator can be added to complete meshes for one unit cell of fuel assembly."
47 : "Elements are quad meshes.");
48 14285 : params.addParam<MooseEnum>("portion", portion, "Control of which part of mesh is created");
49 14285 : params.addRequiredParam<bool>(
50 : "preserve_volumes", "Volume of concentric circles can be preserved using this function.");
51 14285 : params.addClassDescription("This ConcentricCircleMesh source code is to generate concentric "
52 : "circle meshes.");
53 28570 : return params;
54 14285 : }
55 :
56 10 : ConcentricCircleMesh::ConcentricCircleMesh(const InputParameters & parameters)
57 : : MooseMesh(parameters),
58 10 : _num_sectors(getParam<unsigned int>("num_sectors")),
59 10 : _radii(getParam<std::vector<Real>>("radii")),
60 10 : _rings(getParam<std::vector<unsigned int>>("rings")),
61 10 : _inner_mesh_fraction(getParam<Real>("inner_mesh_fraction")),
62 10 : _has_outer_square(getParam<bool>("has_outer_square")),
63 10 : _pitch(getParam<Real>("pitch")),
64 10 : _preserve_volumes(getParam<bool>("preserve_volumes")),
65 20 : _portion(getParam<MooseEnum>("portion"))
66 : {
67 :
68 10 : if (_num_sectors % 2 != 0)
69 0 : mooseError("ConcentricCircleMesh: num_sectors must be an even number.");
70 :
71 : // radii data check
72 80 : for (unsigned i = 0; i < _radii.size() - 1; ++i)
73 70 : if (_radii[i] > _radii[i + 1])
74 0 : mooseError("Radii must be provided in order by starting with the smallest radius and "
75 : "providing the following gradual radii.");
76 :
77 : // condition for setting the size of inner squares.
78 10 : if (_inner_mesh_fraction > std::cos(M_PI / 4))
79 0 : mooseError("The aspect ratio can not be larger than cos(PI/4).");
80 :
81 : // size of 'rings' check
82 10 : if (_has_outer_square)
83 : {
84 10 : if (_rings.size() != _radii.size() + 1)
85 0 : mooseError("The size of 'rings' must be equal to the size of 'radii' plus 1.");
86 : }
87 : else
88 : {
89 0 : if (_rings.size() != _radii.size())
90 0 : mooseError("The size of 'rings' must be equal to the size of 'radii'.");
91 : }
92 : // pitch / 2 must be bigger than any raddi.
93 10 : if (_has_outer_square)
94 90 : for (unsigned i = 0; i < _radii.size(); ++i)
95 80 : if (_pitch / 2 < _radii[i])
96 0 : mooseError("The pitch / 2 must be larger than any radii.");
97 10 : }
98 :
99 : std::unique_ptr<MooseMesh>
100 0 : ConcentricCircleMesh::safeClone() const
101 : {
102 0 : return _app.getFactory().copyConstruct(*this);
103 : }
104 :
105 : void
106 9 : ConcentricCircleMesh::buildMesh()
107 : {
108 : // Get the actual libMesh mesh
109 9 : ReplicatedMesh & mesh = cast_ref<ReplicatedMesh &>(getMesh());
110 : // Set dimension of mesh
111 9 : mesh.set_mesh_dimension(2);
112 9 : mesh.set_spatial_dimension(2);
113 9 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
114 :
115 : // Creating real mesh concentric circles
116 : // i: index for _rings, j: index for _radii
117 9 : std::vector<Real> total_concentric_circles;
118 9 : unsigned int j = 0;
119 81 : while (j < _radii.size())
120 : {
121 72 : unsigned int i = 0;
122 72 : if (j == 0)
123 99 : while (i < _rings[j])
124 : {
125 90 : total_concentric_circles.push_back(_inner_mesh_fraction * _radii[j] +
126 90 : (_radii[j] - _inner_mesh_fraction * _radii[j]) /
127 90 : _rings[j] * (i + 1));
128 90 : ++i;
129 : }
130 : else
131 315 : while (i < _rings[j])
132 : {
133 252 : total_concentric_circles.push_back(_radii[j - 1] +
134 252 : (_radii[j] - _radii[j - 1]) / _rings[j] * (i + 1));
135 252 : ++i;
136 : }
137 72 : ++j;
138 : }
139 :
140 : // volume preserving function is used to conserve volume.
141 9 : const Real d_angle = M_PI / 2 / _num_sectors;
142 :
143 9 : if (_preserve_volumes)
144 : {
145 0 : Real original_radius = 0.0;
146 0 : for (unsigned i = 0; i < total_concentric_circles.size(); ++i)
147 : {
148 : // volume preserving function for the center circle
149 0 : if (i == 0)
150 : {
151 0 : const Real target_area = M_PI * libMesh::Utility::pow<2>(total_concentric_circles[i]);
152 0 : Real modified_radius = std::sqrt(2 * target_area / std::sin(d_angle) / _num_sectors / 4);
153 0 : original_radius = total_concentric_circles[i];
154 0 : total_concentric_circles[i] = modified_radius;
155 : }
156 : else
157 : {
158 : // volume preserving functions for outer circles
159 0 : const Real target_area = M_PI * (libMesh::Utility::pow<2>(total_concentric_circles[i]) -
160 0 : libMesh::Utility::pow<2>(original_radius));
161 0 : Real modified_radius = std::sqrt(target_area / std::sin(d_angle) / _num_sectors / 2 +
162 0 : libMesh::Utility::pow<2>(total_concentric_circles[i - 1]));
163 0 : original_radius = total_concentric_circles[i];
164 0 : total_concentric_circles[i] = modified_radius;
165 : }
166 : }
167 : }
168 :
169 : // number of total nodes
170 9 : unsigned num_total_nodes = 0;
171 9 : if (_has_outer_square)
172 9 : num_total_nodes = libMesh::Utility::pow<2>(_num_sectors / 2 + 1) +
173 9 : (_num_sectors + 1) * (total_concentric_circles.size() + _rings.back()) +
174 9 : (_num_sectors + 1);
175 : else
176 0 : num_total_nodes = libMesh::Utility::pow<2>(_num_sectors / 2 + 1) +
177 0 : (_num_sectors + 1) * total_concentric_circles.size();
178 :
179 9 : std::vector<Node *> nodes(num_total_nodes);
180 9 : unsigned node_id = 0;
181 :
182 : // for adding nodes for the square at the center of the circle
183 45 : for (unsigned i = 0; i <= _num_sectors / 2; ++i)
184 : {
185 36 : const Real x = i * _inner_mesh_fraction * total_concentric_circles[0] / (_num_sectors / 2);
186 180 : for (unsigned j = 0; j <= _num_sectors / 2; ++j)
187 : {
188 144 : const Real y = j * _inner_mesh_fraction * total_concentric_circles[0] / (_num_sectors / 2);
189 144 : nodes[node_id] = mesh.add_point(Point(x, y, 0.0), node_id);
190 144 : ++node_id;
191 : }
192 : }
193 :
194 : // for adding the outer nodes of the square
195 9 : Real current_radius = 0.0;
196 :
197 351 : for (unsigned layers = 0; layers < total_concentric_circles.size(); ++layers)
198 : {
199 342 : current_radius = total_concentric_circles[layers];
200 2736 : for (unsigned num_outer_nodes = 0; num_outer_nodes <= _num_sectors; ++num_outer_nodes)
201 : {
202 2394 : const Real x = current_radius * std::cos(num_outer_nodes * d_angle);
203 2394 : const Real y = current_radius * std::sin(num_outer_nodes * d_angle);
204 2394 : nodes[node_id] = mesh.add_point(Point(x, y, 0.0), node_id);
205 2394 : ++node_id;
206 : }
207 : }
208 :
209 : // adding nodes for the unit cell of fuel assembly.
210 9 : if (_has_outer_square)
211 : {
212 9 : Real current_radius_moderator = 0.0;
213 99 : for (unsigned i = 1; i <= _rings.back(); ++i)
214 : {
215 90 : current_radius_moderator =
216 90 : _radii.back() + i * (_pitch / 2 - _radii.back()) / (_rings.back() + 1);
217 90 : total_concentric_circles.push_back(current_radius_moderator);
218 720 : for (unsigned num_outer_nodes = 0; num_outer_nodes <= _num_sectors; ++num_outer_nodes)
219 : {
220 630 : const Real x = current_radius_moderator * std::cos(num_outer_nodes * d_angle);
221 630 : const Real y = current_radius_moderator * std::sin(num_outer_nodes * d_angle);
222 630 : nodes[node_id] = mesh.add_point(Point(x, y, 0.0), node_id);
223 630 : ++node_id;
224 : }
225 : }
226 :
227 45 : for (unsigned j = 0; j < _num_sectors / 2 + 1; ++j)
228 : {
229 36 : const Real x = _pitch / 2;
230 36 : const Real y = _pitch / 2 * std::tan(j * d_angle);
231 36 : nodes[node_id] = mesh.add_point(Point(x, y, 0.0), node_id);
232 36 : ++node_id;
233 : }
234 :
235 36 : for (unsigned i = 0; i < _num_sectors / 2; ++i)
236 : {
237 27 : const Real x = _pitch / 2 * std::cos((i + _num_sectors / 2 + 1) * d_angle) /
238 27 : std::sin((i + _num_sectors / 2 + 1) * d_angle);
239 27 : const Real y = _pitch / 2;
240 27 : nodes[node_id] = mesh.add_point(Point(x, y, 0.0), node_id);
241 27 : ++node_id;
242 : }
243 : }
244 : // Currently, index, limit, counter variables use the int type because of the 'modulo' function.
245 : // adding elements
246 9 : int index = 0;
247 9 : int limit = 0;
248 9 : int standard = static_cast<int>(_num_sectors);
249 :
250 : // This is to set the limit for the index
251 9 : if (standard > 4)
252 : {
253 9 : int additional_term = 0;
254 9 : int counter = standard;
255 18 : while (counter > 4)
256 : {
257 9 : counter = counter - 2;
258 9 : additional_term = additional_term + counter;
259 : }
260 9 : limit = standard + additional_term;
261 : }
262 0 : else if (standard == 4)
263 0 : limit = standard;
264 :
265 : // SubdomainIDs set up
266 9 : std::vector<unsigned int> subdomainIDs;
267 90 : for (unsigned int i = 0; i < _rings.size(); ++i)
268 513 : for (unsigned int j = 0; j < _rings[i]; ++j)
269 432 : subdomainIDs.push_back(i + 1);
270 :
271 9 : if (_has_outer_square)
272 9 : subdomainIDs.push_back(subdomainIDs.back());
273 : // adding elements in the square
274 90 : while (index <= limit)
275 : {
276 81 : Elem * elem = mesh.add_elem(new Quad4);
277 81 : elem->set_node(0, nodes[index]);
278 81 : elem->set_node(1, nodes[index + _num_sectors / 2 + 1]);
279 81 : elem->set_node(2, nodes[index + _num_sectors / 2 + 2]);
280 81 : elem->set_node(3, nodes[index + 1]);
281 81 : elem->subdomain_id() = subdomainIDs[0];
282 :
283 81 : if (index < standard / 2)
284 27 : boundary_info.add_side(elem, 3, 1);
285 81 : if (index % (standard / 2 + 1) == 0)
286 27 : boundary_info.add_side(elem, 0, 2);
287 :
288 81 : ++index;
289 81 : if ((index - standard / 2) % (standard / 2 + 1) == 0)
290 27 : ++index;
291 : }
292 :
293 9 : index = (_num_sectors / 2 + 1) * (_num_sectors / 2);
294 9 : limit = (_num_sectors / 2) * (_num_sectors / 2 + 2);
295 :
296 : // adding elements in one outer layer of the square (right side)
297 36 : while (index < limit)
298 : {
299 27 : Elem * elem = mesh.add_elem(new Quad4);
300 27 : elem->set_node(0, nodes[index]);
301 27 : elem->set_node(1, nodes[index + _num_sectors / 2 + 1]);
302 27 : elem->set_node(2, nodes[index + _num_sectors / 2 + 2]);
303 27 : elem->set_node(3, nodes[index + 1]);
304 27 : elem->subdomain_id() = subdomainIDs[0];
305 :
306 27 : if (index == (standard / 2 + 1) * (standard / 2))
307 9 : boundary_info.add_side(elem, 0, 2);
308 :
309 27 : ++index;
310 : }
311 :
312 : // adding elements in one outer layer of the square (left side)
313 9 : int counter = 0;
314 36 : while (index != standard / 2)
315 : {
316 27 : Elem * elem = mesh.add_elem(new Quad4);
317 27 : elem->set_node(0, nodes[index]);
318 27 : elem->set_node(1, nodes[index + (_num_sectors / 2 + 1) + counter * (_num_sectors / 2 + 2)]);
319 27 : elem->set_node(2, nodes[index + (_num_sectors / 2 + 1) + counter * (_num_sectors / 2 + 2) + 1]);
320 27 : elem->set_node(3, nodes[index - _num_sectors / 2 - 1]);
321 27 : elem->subdomain_id() = subdomainIDs[0];
322 :
323 27 : if (index == standard + 1)
324 9 : boundary_info.add_side(elem, 2, 1);
325 :
326 27 : index = index - _num_sectors / 2 - 1;
327 27 : ++counter;
328 : }
329 :
330 : // adding elements for other concentric circles
331 9 : index = libMesh::Utility::pow<2>(_num_sectors / 2 + 1);
332 9 : limit = static_cast<int>(num_total_nodes) - standard - 2;
333 9 : int num_nodes_boundary = libMesh::Utility::pow<2>(_num_sectors / 2 + 1) + _num_sectors + 1;
334 :
335 9 : counter = 0;
336 2601 : while (index < limit)
337 : {
338 :
339 2592 : Elem * elem = mesh.add_elem(new Quad4);
340 2592 : elem->set_node(0, nodes[index]);
341 2592 : elem->set_node(1, nodes[index + _num_sectors + 1]);
342 2592 : elem->set_node(2, nodes[index + _num_sectors + 2]);
343 2592 : elem->set_node(3, nodes[index + 1]);
344 :
345 127008 : for (int i = 0; i < static_cast<int>(subdomainIDs.size()) - 1; ++i)
346 124416 : if (index < limit - (standard + 1) * i && index >= limit - (standard + 1) * (i + 1))
347 2592 : elem->subdomain_id() = subdomainIDs[subdomainIDs.size() - 1 - i];
348 :
349 2592 : int const initial = libMesh::Utility::pow<2>(standard / 2 + 1);
350 2592 : int const final = libMesh::Utility::pow<2>(standard / 2 + 1) + _num_sectors - 1;
351 :
352 2592 : if ((index - initial) % (standard + 1) == 0)
353 432 : boundary_info.add_side(elem, 0, 2);
354 2592 : if ((index - final) % (standard + 1) == 0)
355 432 : boundary_info.add_side(elem, 2, 1);
356 2592 : if (index > limit - (standard + 1))
357 : {
358 54 : if (_has_outer_square)
359 : {
360 54 : if (index < limit - standard + standard / 2)
361 27 : boundary_info.add_side(elem, 1, 3);
362 : else
363 27 : boundary_info.add_side(elem, 1, 4);
364 : }
365 : else
366 : {
367 0 : boundary_info.add_side(elem, 1, 3);
368 : }
369 : }
370 2592 : ++index;
371 2592 : if (index == (num_nodes_boundary + counter * (standard + 1)) - 1)
372 : {
373 432 : ++index;
374 432 : ++counter;
375 : }
376 : }
377 :
378 : // This is to set boundary names.
379 9 : boundary_info.sideset_name(1) = "left";
380 9 : boundary_info.sideset_name(2) = "bottom";
381 :
382 9 : if (!_has_outer_square)
383 0 : boundary_info.sideset_name(3) = "outer";
384 : else
385 : {
386 9 : boundary_info.sideset_name(3) = "right";
387 9 : boundary_info.sideset_name(4) = "top";
388 : }
389 :
390 9 : if (_portion == "top_left")
391 : {
392 0 : libMesh::MeshTools::Modification::rotate(mesh, 90, 0, 0);
393 0 : boundary_info.sideset_name(1) = "bottom";
394 0 : boundary_info.sideset_name(2) = "right";
395 :
396 0 : if (!_has_outer_square)
397 0 : boundary_info.sideset_name(3) = "outer";
398 : else
399 : {
400 0 : boundary_info.sideset_name(3) = "top";
401 0 : boundary_info.sideset_name(4) = "left";
402 : }
403 : }
404 9 : else if (_portion == "bottom_left")
405 : {
406 0 : libMesh::MeshTools::Modification::rotate(mesh, 180, 0, 0);
407 0 : boundary_info.sideset_name(1) = "right";
408 0 : boundary_info.sideset_name(2) = "top";
409 :
410 0 : if (!_has_outer_square)
411 0 : boundary_info.sideset_name(3) = "outer";
412 : else
413 : {
414 0 : boundary_info.sideset_name(3) = "left";
415 0 : boundary_info.sideset_name(4) = "bottom";
416 : }
417 : }
418 9 : else if (_portion == "bottom_right")
419 : {
420 0 : libMesh::MeshTools::Modification::rotate(mesh, 270, 0, 0);
421 0 : boundary_info.sideset_name(1) = "top";
422 0 : boundary_info.sideset_name(2) = "left";
423 :
424 0 : if (!_has_outer_square)
425 0 : boundary_info.sideset_name(3) = "outer";
426 : else
427 : {
428 0 : boundary_info.sideset_name(3) = "bottom";
429 0 : boundary_info.sideset_name(4) = "right";
430 : }
431 : }
432 :
433 9 : else if (_portion == "top_half")
434 : {
435 0 : ReplicatedMesh other_mesh(mesh);
436 : // This is to rotate the mesh and also to reset boundary IDs.
437 0 : libMesh::MeshTools::Modification::rotate(other_mesh, 90, 0, 0);
438 0 : if (_has_outer_square)
439 : {
440 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
441 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 2, 6);
442 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 3, 7);
443 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 4, 1);
444 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 5, 2);
445 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 6, 3);
446 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 7, 4);
447 0 : mesh.prepare_for_use();
448 0 : other_mesh.prepare_for_use();
449 0 : mesh.stitch_meshes(other_mesh, 1, 3, TOLERANCE, true);
450 0 : mesh.get_boundary_info().sideset_name(1) = "left";
451 0 : mesh.get_boundary_info().sideset_name(2) = "bottom";
452 0 : mesh.get_boundary_info().sideset_name(3) = "right";
453 0 : mesh.get_boundary_info().sideset_name(4) = "top";
454 : }
455 : else
456 : {
457 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
458 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 2, 1);
459 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 5, 2);
460 0 : mesh.prepare_for_use();
461 0 : other_mesh.prepare_for_use();
462 0 : mesh.stitch_meshes(other_mesh, 1, 1, TOLERANCE, true);
463 :
464 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 2, 1);
465 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 3, 2);
466 0 : mesh.get_boundary_info().sideset_name(1) = "bottom";
467 0 : mesh.get_boundary_info().sideset_name(2) = "outer";
468 : }
469 0 : other_mesh.clear();
470 0 : }
471 :
472 9 : else if (_portion == "right_half")
473 : {
474 0 : ReplicatedMesh other_mesh(mesh);
475 : // This is to rotate the mesh and also to reset boundary IDs.
476 0 : libMesh::MeshTools::Modification::rotate(other_mesh, 270, 0, 0);
477 0 : if (_has_outer_square)
478 : {
479 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
480 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 2, 6);
481 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 3, 7);
482 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 4, 3);
483 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 5, 4);
484 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 6, 1);
485 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 7, 2);
486 0 : mesh.prepare_for_use();
487 0 : other_mesh.prepare_for_use();
488 0 : mesh.stitch_meshes(other_mesh, 2, 4, TOLERANCE, true);
489 0 : mesh.get_boundary_info().sideset_name(1) = "left";
490 0 : mesh.get_boundary_info().sideset_name(2) = "bottom";
491 0 : mesh.get_boundary_info().sideset_name(3) = "right";
492 0 : mesh.get_boundary_info().sideset_name(4) = "top";
493 : }
494 : else
495 : {
496 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
497 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 2, 1);
498 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 5, 2);
499 0 : mesh.prepare_for_use();
500 0 : other_mesh.prepare_for_use();
501 0 : mesh.stitch_meshes(other_mesh, 2, 2, TOLERANCE, true);
502 :
503 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 3, 2);
504 0 : mesh.get_boundary_info().sideset_name(1) = "left";
505 0 : mesh.get_boundary_info().sideset_name(2) = "outer";
506 : }
507 0 : other_mesh.clear();
508 0 : }
509 9 : else if (_portion == "left_half")
510 : {
511 0 : ReplicatedMesh other_mesh(mesh);
512 :
513 : // This is to rotate the mesh and to reset boundary IDs.
514 0 : libMesh::MeshTools::Modification::rotate(other_mesh, 90, 0, 0);
515 0 : libMesh::MeshTools::Modification::rotate(mesh, 180, 0, 0);
516 0 : if (_has_outer_square)
517 : {
518 : // The other mesh is created by rotating the original mesh about 90 degrees.
519 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
520 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 2, 6);
521 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 3, 7);
522 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 4, 1);
523 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 5, 2);
524 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 6, 3);
525 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 7, 4);
526 : // The original mesh is then rotated about 180 degrees.
527 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 1, 5);
528 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 2, 6);
529 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 3, 7);
530 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 4, 2);
531 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 5, 3);
532 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 6, 4);
533 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 7, 1);
534 0 : mesh.prepare_for_use();
535 0 : other_mesh.prepare_for_use();
536 0 : mesh.stitch_meshes(other_mesh, 4, 2, TOLERANCE, true);
537 0 : mesh.get_boundary_info().sideset_name(1) = "left";
538 0 : mesh.get_boundary_info().sideset_name(2) = "bottom";
539 0 : mesh.get_boundary_info().sideset_name(3) = "right";
540 0 : mesh.get_boundary_info().sideset_name(4) = "top";
541 : }
542 : else
543 : {
544 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 1, 5);
545 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 2, 1);
546 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 5, 2);
547 0 : mesh.prepare_for_use();
548 0 : other_mesh.prepare_for_use();
549 0 : mesh.stitch_meshes(other_mesh, 1, 1, TOLERANCE, true);
550 :
551 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 2, 1);
552 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 3, 2);
553 0 : mesh.get_boundary_info().sideset_name(1) = "right";
554 0 : mesh.get_boundary_info().sideset_name(2) = "outer";
555 : }
556 0 : other_mesh.clear();
557 0 : }
558 9 : else if (_portion == "bottom_half")
559 : {
560 0 : ReplicatedMesh other_mesh(mesh);
561 : // This is to rotate the mesh and also to reset boundary IDs.
562 0 : libMesh::MeshTools::Modification::rotate(other_mesh, 180, 0, 0);
563 0 : libMesh::MeshTools::Modification::rotate(mesh, 270, 0, 0);
564 0 : if (_has_outer_square)
565 : {
566 : // The other mesh is created by rotating the original mesh about 180 degrees.
567 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
568 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 2, 6);
569 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 3, 7);
570 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 4, 2);
571 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 5, 3);
572 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 6, 4);
573 0 : libMesh::MeshTools::Modification::change_boundary_id(other_mesh, 7, 1);
574 : // The original mesh is rotated about 270 degrees.
575 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 1, 5);
576 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 2, 6);
577 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 3, 7);
578 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 4, 3);
579 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 5, 4);
580 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 6, 1);
581 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 7, 2);
582 0 : mesh.prepare_for_use();
583 0 : other_mesh.prepare_for_use();
584 0 : mesh.stitch_meshes(other_mesh, 1, 3, TOLERANCE, true);
585 0 : mesh.get_boundary_info().sideset_name(1) = "left";
586 0 : mesh.get_boundary_info().sideset_name(2) = "bottom";
587 0 : mesh.get_boundary_info().sideset_name(3) = "right";
588 0 : mesh.get_boundary_info().sideset_name(4) = "top";
589 : }
590 : else
591 : {
592 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 1, 5);
593 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 2, 1);
594 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 5, 2);
595 0 : mesh.prepare_for_use();
596 0 : other_mesh.prepare_for_use();
597 0 : mesh.stitch_meshes(other_mesh, 1, 1, TOLERANCE, true);
598 :
599 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 2, 1);
600 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 3, 2);
601 0 : mesh.get_boundary_info().sideset_name(1) = "top";
602 0 : mesh.get_boundary_info().sideset_name(2) = "outer";
603 : }
604 0 : other_mesh.clear();
605 0 : }
606 9 : else if (_portion == "full")
607 : {
608 9 : ReplicatedMesh portion_two(mesh);
609 :
610 : // This is to rotate the mesh and also to reset boundary IDs.
611 9 : libMesh::MeshTools::Modification::rotate(portion_two, 90, 0, 0);
612 :
613 9 : if (_has_outer_square)
614 : {
615 : // Portion 2: 2nd quadrant
616 9 : libMesh::MeshTools::Modification::change_boundary_id(portion_two, 1, 5);
617 9 : libMesh::MeshTools::Modification::change_boundary_id(portion_two, 2, 6);
618 9 : libMesh::MeshTools::Modification::change_boundary_id(portion_two, 3, 7);
619 9 : libMesh::MeshTools::Modification::change_boundary_id(portion_two, 4, 1);
620 9 : libMesh::MeshTools::Modification::change_boundary_id(portion_two, 5, 2);
621 9 : libMesh::MeshTools::Modification::change_boundary_id(portion_two, 6, 3);
622 9 : libMesh::MeshTools::Modification::change_boundary_id(portion_two, 7, 4);
623 9 : mesh.prepare_for_use();
624 9 : portion_two.prepare_for_use();
625 : // 'top_half'
626 9 : mesh.stitch_meshes(portion_two, 1, 3, TOLERANCE, true);
627 :
628 : // 'bottom_half'
629 9 : ReplicatedMesh portion_bottom(mesh);
630 9 : libMesh::MeshTools::Modification::rotate(portion_bottom, 180, 0, 0);
631 9 : libMesh::MeshTools::Modification::change_boundary_id(portion_bottom, 1, 5);
632 9 : libMesh::MeshTools::Modification::change_boundary_id(portion_bottom, 2, 6);
633 9 : libMesh::MeshTools::Modification::change_boundary_id(portion_bottom, 3, 7);
634 9 : libMesh::MeshTools::Modification::change_boundary_id(portion_bottom, 4, 2);
635 9 : libMesh::MeshTools::Modification::change_boundary_id(portion_bottom, 5, 3);
636 9 : libMesh::MeshTools::Modification::change_boundary_id(portion_bottom, 6, 4);
637 9 : libMesh::MeshTools::Modification::change_boundary_id(portion_bottom, 7, 1);
638 9 : mesh.prepare_for_use();
639 9 : portion_bottom.prepare_for_use();
640 : // 'full'
641 9 : mesh.stitch_meshes(portion_bottom, 2, 4, TOLERANCE, true);
642 :
643 9 : mesh.get_boundary_info().sideset_name(1) = "left";
644 9 : mesh.get_boundary_info().sideset_name(2) = "bottom";
645 9 : mesh.get_boundary_info().sideset_name(3) = "right";
646 9 : mesh.get_boundary_info().sideset_name(4) = "top";
647 9 : portion_bottom.clear();
648 9 : }
649 : else
650 : {
651 0 : libMesh::MeshTools::Modification::change_boundary_id(portion_two, 1, 5);
652 0 : libMesh::MeshTools::Modification::change_boundary_id(portion_two, 2, 1);
653 0 : libMesh::MeshTools::Modification::change_boundary_id(portion_two, 5, 2);
654 : // 'top half'
655 0 : mesh.prepare_for_use();
656 0 : portion_two.prepare_for_use();
657 0 : mesh.stitch_meshes(portion_two, 1, 1, TOLERANCE, true);
658 : // 'bottom half'
659 0 : ReplicatedMesh portion_bottom(mesh);
660 0 : libMesh::MeshTools::Modification::rotate(portion_bottom, 180, 0, 0);
661 : // 'full'
662 0 : mesh.prepare_for_use();
663 0 : portion_bottom.prepare_for_use();
664 0 : mesh.stitch_meshes(portion_bottom, 2, 2, TOLERANCE, true);
665 0 : libMesh::MeshTools::Modification::change_boundary_id(mesh, 3, 1);
666 0 : mesh.get_boundary_info().sideset_name(1) = "outer";
667 0 : portion_bottom.clear();
668 0 : }
669 9 : portion_two.clear();
670 9 : }
671 18 : if (_portion != "top_half" && _portion != "right_half" && _portion != "left_half" &&
672 18 : _portion != "bottom_half" && _portion != "full")
673 0 : mesh.prepare_for_use();
674 9 : }
|
