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 "ConcentricCircleMeshGenerator.h"
11 : #include "CastUniquePointer.h"
12 : #include "libmesh/replicated_mesh.h"
13 : #include "libmesh/face_quad4.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 : #include "libmesh/mesh_smoother_laplace.h"
19 :
20 : // C++ includes
21 : #include <cmath> // provides round, not std::round (see http://www.cplusplus.com/reference/cmath/round/)
22 :
23 : registerMooseObject("MooseApp", ConcentricCircleMeshGenerator);
24 :
25 : InputParameters
26 3839 : ConcentricCircleMeshGenerator::validParams()
27 : {
28 3839 : InputParameters params = MeshGenerator::validParams();
29 15356 : params.addParam<SubdomainName>("subdomain_name",
30 : "Name of the subdomain assigned to all the elements");
31 : MooseEnum portion(
32 : "full top_right top_left bottom_left bottom_right right_half left_half top_half bottom_half",
33 15356 : "full");
34 15356 : params.addRequiredParam<unsigned int>("num_sectors",
35 : "num_sectors % 2 = 0, num_sectors > 0"
36 : "Number of azimuthal sectors in each quadrant"
37 : "'num_sectors' must be an even number.");
38 15356 : params.addRequiredParam<std::vector<Real>>("radii", "Radii of major concentric circles");
39 15356 : params.addRequiredParam<std::vector<unsigned int>>(
40 : "rings", "Number of rings in each circle or in the enclosing square");
41 15356 : params.addRequiredParam<bool>(
42 : "has_outer_square",
43 : "It determines if meshes for a outer square are added to concentric circle meshes.");
44 19195 : params.addRangeCheckedParam<Real>(
45 : "pitch",
46 7678 : 0.0,
47 : "pitch>=0.0",
48 : "The enclosing square can be added to the completed concentric circle mesh."
49 : "Elements are quad meshes.");
50 15356 : params.addParam<MooseEnum>("portion", portion, "Control of which part of mesh is created");
51 15356 : params.addRequiredParam<bool>(
52 : "preserve_volumes", "Volume of concentric circles can be preserved using this function.");
53 15356 : params.addParam<unsigned int>("smoothing_max_it", 1, "Number of Laplacian smoothing iterations");
54 3839 : params.addClassDescription(
55 : "This ConcentricCircleMeshGenerator source code is to generate concentric "
56 : "circle meshes.");
57 :
58 7678 : return params;
59 3839 : }
60 :
61 371 : ConcentricCircleMeshGenerator::ConcentricCircleMeshGenerator(const InputParameters & parameters)
62 : : MeshGenerator(parameters),
63 371 : _num_sectors(getParam<unsigned int>("num_sectors")),
64 742 : _radii(getParam<std::vector<Real>>("radii")),
65 742 : _rings(getParam<std::vector<unsigned int>>("rings")),
66 742 : _has_outer_square(getParam<bool>("has_outer_square")),
67 742 : _pitch(getParam<Real>("pitch")),
68 742 : _preserve_volumes(getParam<bool>("preserve_volumes")),
69 742 : _smoothing_max_it(getParam<unsigned int>("smoothing_max_it")),
70 1113 : _portion(getParam<MooseEnum>("portion"))
71 : {
72 742 : declareMeshProperty("use_distributed_mesh", false);
73 :
74 371 : if (_num_sectors % 2 != 0)
75 0 : paramError("num_sectors", "must be an even number.");
76 :
77 : // radii data check
78 723 : for (unsigned i = 0; i < _radii.size() - 1; ++i)
79 : {
80 352 : if (_radii[i] == 0.0)
81 0 : paramError("radii", "must be positive numbers.");
82 352 : if (_radii[i] > _radii[i + 1])
83 0 : paramError("radii",
84 : "must be provided in order by starting with the smallest radius providing the "
85 : "following gradual radii.");
86 : }
87 :
88 : // size of 'rings' check
89 371 : if (_has_outer_square)
90 : {
91 104 : if (_rings.size() != _radii.size() + 1)
92 0 : paramError("rings", "The size of 'rings' must be equal to the size of 'radii' plus 1.");
93 : }
94 : else
95 : {
96 267 : if (_rings.size() != _radii.size())
97 0 : paramError("rings", "The size of 'rings' must be equal to the size of 'radii'.");
98 : }
99 : // pitch / 2 must be bigger than any raddi.
100 371 : if (_has_outer_square)
101 444 : for (unsigned i = 0; i < _radii.size(); ++i)
102 340 : if (_pitch / 2 < _radii[i])
103 0 : paramError("pitch", "The pitch / 2 must be larger than any radii.");
104 371 : }
105 :
106 : std::unique_ptr<MeshBase>
107 369 : ConcentricCircleMeshGenerator::generate()
108 : {
109 369 : auto mesh = buildReplicatedMesh(2);
110 369 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
111 :
112 : // Creating real mesh concentric circles
113 : // i: index for _rings, j: index for _radii
114 369 : std::vector<Real> total_concentric_circles;
115 369 : unsigned int j = 0;
116 1082 : while (j < _radii.size())
117 : {
118 713 : unsigned int i = 0;
119 713 : if (j == 0)
120 1121 : while (i < _rings[j])
121 : {
122 752 : total_concentric_circles.push_back(_radii[j] / (_num_sectors / 2 + _rings[j]) *
123 752 : (i + _num_sectors / 2 + 1));
124 752 : ++i;
125 : }
126 : else
127 1038 : while (i < _rings[j])
128 : {
129 694 : total_concentric_circles.push_back(_radii[j - 1] +
130 694 : (_radii[j] - _radii[j - 1]) / _rings[j] * (i + 1));
131 694 : ++i;
132 : }
133 713 : ++j;
134 : }
135 :
136 : // volume preserving function is used to conserve volume.
137 369 : const Real d_angle = M_PI / 2 / _num_sectors;
138 :
139 369 : if (_preserve_volumes)
140 : {
141 138 : Real original_radius = 0.0;
142 564 : for (unsigned i = 0; i < total_concentric_circles.size(); ++i)
143 : {
144 : // volume preserving function for the center circle
145 426 : if (i == 0)
146 : {
147 138 : const Real target_area = M_PI * Utility::pow<2>(total_concentric_circles[i]);
148 138 : Real modified_radius = std::sqrt(2 * target_area / std::sin(d_angle) / _num_sectors / 4);
149 138 : original_radius = total_concentric_circles[i];
150 138 : total_concentric_circles[i] = modified_radius;
151 : }
152 : else
153 : {
154 : // volume preserving functions for outer circles
155 288 : const Real target_area = M_PI * (Utility::pow<2>(total_concentric_circles[i]) -
156 288 : Utility::pow<2>(original_radius));
157 576 : Real modified_radius = std::sqrt(target_area / std::sin(d_angle) / _num_sectors / 2 +
158 288 : Utility::pow<2>(total_concentric_circles[i - 1]));
159 288 : original_radius = total_concentric_circles[i];
160 288 : total_concentric_circles[i] = modified_radius;
161 : }
162 : }
163 : }
164 :
165 : // number of total nodes
166 369 : unsigned num_total_nodes = 0;
167 :
168 369 : if (_has_outer_square)
169 104 : num_total_nodes = Utility::pow<2>(_num_sectors / 2 + 1) +
170 104 : (_num_sectors + 1) * total_concentric_circles.size() +
171 104 : Utility::pow<2>(_rings.back() + 2) + _num_sectors * (_rings.back() + 1) - 1;
172 : else
173 530 : num_total_nodes = Utility::pow<2>(_num_sectors / 2 + 1) +
174 265 : (_num_sectors + 1) * total_concentric_circles.size();
175 :
176 369 : std::vector<Node *> nodes(num_total_nodes);
177 :
178 369 : unsigned node_id = 0;
179 :
180 : // for adding nodes for the square at the center of the circle
181 369 : Real xx = total_concentric_circles[0] / (_num_sectors / 2 + 1) * _num_sectors / 2;
182 369 : Point p1 = Point(xx, 0, 0);
183 369 : Point p2 = Point(0, xx, 0);
184 369 : Point p3 = Point(xx * std::sqrt(2.0) / 2, xx * std::sqrt(2.0) / 2, 0);
185 1781 : for (unsigned i = 0; i <= _num_sectors / 2; ++i)
186 : {
187 1412 : Real fx = i / (_num_sectors / 2.0);
188 7164 : for (unsigned j = 0; j <= _num_sectors / 2; ++j)
189 : {
190 5752 : Real fy = j / (_num_sectors / 2.0);
191 5752 : Point p = p1 * fx * (1 - fy) + p2 * fy * (1 - fx) + p3 * fx * fy;
192 5752 : nodes[node_id] = mesh->add_point(p, node_id);
193 5752 : ++node_id;
194 : }
195 : }
196 :
197 : // for adding the outer layers of the square
198 369 : Real current_radius = total_concentric_circles[0];
199 :
200 : // for adding the outer circles of the square.
201 1815 : for (unsigned layers = 0; layers < total_concentric_circles.size(); ++layers)
202 : {
203 1446 : current_radius = total_concentric_circles[layers];
204 11756 : for (unsigned num_outer_nodes = 0; num_outer_nodes <= _num_sectors; ++num_outer_nodes)
205 : {
206 10310 : const Real x = current_radius * std::cos(num_outer_nodes * d_angle);
207 10310 : const Real y = current_radius * std::sin(num_outer_nodes * d_angle);
208 10310 : nodes[node_id] = mesh->add_point(Point(x, y, 0.0), node_id);
209 10310 : ++node_id;
210 : }
211 : }
212 :
213 : // adding nodes for the enclosing square.
214 : // adding nodes for the left edge of the square.
215 : // applying the method for partitioning the line segments.
216 :
217 369 : if (_has_outer_square)
218 : {
219 : // putting nodes for the left side of the enclosing square.
220 460 : for (unsigned i = 0; i <= _num_sectors / 2; ++i)
221 : {
222 356 : const Real a1 = (_pitch / 2) * i / (_num_sectors / 2 + _rings.back() + 1);
223 356 : const Real b1 = _pitch / 2;
224 :
225 356 : const Real a2 = total_concentric_circles.back() * std::cos(M_PI / 2 - i * d_angle);
226 356 : const Real b2 = total_concentric_circles.back() * std::sin(M_PI / 2 - i * d_angle);
227 :
228 1476 : for (unsigned j = 0; j <= _rings.back(); ++j)
229 : {
230 1120 : Real x = ((_rings.back() + 1 - j) * a1 + j * a2) / (_rings.back() + 1);
231 1120 : Real y = ((_rings.back() + 1 - j) * b1 + j * b2) / (_rings.back() + 1);
232 1120 : nodes[node_id] = mesh->add_point(Point(x, y, 0.0), node_id);
233 1120 : ++node_id;
234 : }
235 : }
236 : // putting nodes for the center part of the enclosing square.
237 104 : unsigned k = 1;
238 414 : for (unsigned i = 1; i <= _rings.back() + 1; ++i)
239 : {
240 : const Real a1 =
241 310 : (_pitch / 2) * (i + _num_sectors / 2) / (_num_sectors / 2 + _rings.back() + 1);
242 310 : const Real b1 = _pitch / 2;
243 :
244 310 : const Real a2 = _pitch / 2;
245 310 : const Real b2 = (_pitch / 2) * (_num_sectors / 2) / (_num_sectors / 2 + _rings.back() + 1);
246 :
247 310 : const Real a3 = total_concentric_circles.back() * std::cos(M_PI / 4);
248 310 : const Real b3 = total_concentric_circles.back() * std::sin(M_PI / 4);
249 :
250 310 : const Real a4 = ((_rings.back() + 1 - k) * a3 + k * a2) / (_rings.back() + 1);
251 310 : const Real b4 = ((_rings.back() + 1 - k) * b3 + k * b2) / (_rings.back() + 1);
252 :
253 1738 : for (unsigned j = 0; j <= _rings.back() + 1; ++j)
254 : {
255 1428 : Real x = ((_rings.back() + 1 - j) * a1 + j * a4) / (_rings.back() + 1);
256 1428 : Real y = ((_rings.back() + 1 - j) * b1 + j * b4) / (_rings.back() + 1);
257 1428 : nodes[node_id] = mesh->add_point(Point(x, y, 0.0), node_id);
258 1428 : ++node_id;
259 : }
260 310 : ++k;
261 : }
262 :
263 : // putting nodes for the right part of the enclosing square.
264 356 : for (unsigned i = 1; i <= _num_sectors / 2; ++i)
265 : {
266 252 : const Real a1 = _pitch / 2;
267 : const Real b1 =
268 252 : (_pitch / 2) * (_num_sectors / 2 - i) / (_num_sectors / 2 + _rings.back() + 1);
269 :
270 252 : const Real a2 = total_concentric_circles.back() * std::cos(M_PI / 4 - i * d_angle);
271 252 : const Real b2 = total_concentric_circles.back() * std::sin(M_PI / 4 - i * d_angle);
272 :
273 1062 : for (unsigned j = 0; j <= _rings.back(); ++j)
274 : {
275 810 : Real x = ((_rings.back() + 1 - j) * a1 + j * a2) / (_rings.back() + 1);
276 810 : Real y = ((_rings.back() + 1 - j) * b1 + j * b2) / (_rings.back() + 1);
277 810 : nodes[node_id] = mesh->add_point(Point(x, y, 0.0), node_id);
278 810 : ++node_id;
279 : }
280 : }
281 : }
282 :
283 : // Currently, index, limit, counter variables use the int type because of the 'modulo' function.
284 : // adding elements
285 369 : int index = 0;
286 369 : int limit = 0;
287 369 : int standard = static_cast<int>(_num_sectors);
288 :
289 : // This is to set the limit for the index
290 369 : if (standard > 4)
291 : {
292 249 : int additional_term = 0;
293 249 : int counter = standard;
294 600 : while (counter > 4)
295 : {
296 351 : counter = counter - 2;
297 351 : additional_term = additional_term + counter;
298 : }
299 249 : limit = standard + additional_term;
300 : }
301 120 : else if (standard == 4)
302 74 : limit = standard;
303 :
304 : // SubdomainIDs set up
305 369 : std::vector<unsigned int> subdomainIDs;
306 :
307 369 : if (_has_outer_square)
308 444 : for (unsigned int i = 0; i < _rings.size() - 1; ++i)
309 996 : for (unsigned int j = 0; j < _rings[i]; ++j)
310 656 : subdomainIDs.push_back(i + 1);
311 : else
312 638 : for (unsigned int i = 0; i < _rings.size(); ++i)
313 1163 : for (unsigned int j = 0; j < _rings[i]; ++j)
314 790 : subdomainIDs.push_back(i + 1);
315 :
316 : // adding elements in the square
317 3666 : while (index <= limit)
318 : {
319 : // inner circle area (polygonal core)
320 3297 : Elem * elem = mesh->add_elem(std::make_unique<Quad4>());
321 3297 : elem->set_node(0, nodes[index]);
322 3297 : elem->set_node(1, nodes[index + _num_sectors / 2 + 1]);
323 3297 : elem->set_node(2, nodes[index + _num_sectors / 2 + 2]);
324 3297 : elem->set_node(3, nodes[index + 1]);
325 3297 : elem->subdomain_id() = subdomainIDs[0];
326 :
327 3297 : if (index < standard / 2)
328 1043 : boundary_info.add_side(elem, 3, 1);
329 3297 : if (index % (standard / 2 + 1) == 0)
330 1043 : boundary_info.add_side(elem, 0, 2);
331 :
332 3297 : ++index;
333 :
334 : // increment by 2 indices is necessary depending on where the index points to.
335 3297 : if ((index - standard / 2) % (standard / 2 + 1) == 0)
336 1043 : ++index;
337 : }
338 :
339 369 : index = (_num_sectors / 2 + 1) * (_num_sectors / 2);
340 369 : limit = (_num_sectors / 2) * (_num_sectors / 2 + 2);
341 :
342 : // adding elements in one outer layer of the square (right side)
343 1412 : while (index < limit)
344 : {
345 : // inner circle elements touching B
346 1043 : Elem * elem = mesh->add_elem(std::make_unique<Quad4>());
347 1043 : elem->set_node(0, nodes[index]);
348 1043 : elem->set_node(1, nodes[index + _num_sectors / 2 + 1]);
349 1043 : elem->set_node(2, nodes[index + _num_sectors / 2 + 2]);
350 1043 : elem->set_node(3, nodes[index + 1]);
351 1043 : elem->subdomain_id() = subdomainIDs[0];
352 :
353 1043 : if (index == (standard / 2 + 1) * (standard / 2))
354 369 : boundary_info.add_side(elem, 0, 2);
355 :
356 1043 : ++index;
357 : }
358 :
359 : // adding elements in one outer layer of the square (left side)
360 369 : int counter = 0;
361 1412 : while (index != standard / 2)
362 : {
363 : // inner circle elements touching C
364 1043 : Elem * elem = mesh->add_elem(std::make_unique<Quad4>());
365 1043 : elem->set_node(0, nodes[index]);
366 1043 : elem->set_node(1, nodes[index + (_num_sectors / 2 + 1) + counter * (_num_sectors / 2 + 2)]);
367 1043 : elem->set_node(2, nodes[index + (_num_sectors / 2 + 1) + counter * (_num_sectors / 2 + 2) + 1]);
368 1043 : elem->set_node(3, nodes[index - _num_sectors / 2 - 1]);
369 1043 : elem->subdomain_id() = subdomainIDs[0];
370 :
371 1043 : if (index == standard + 1)
372 369 : boundary_info.add_side(elem, 2, 1);
373 :
374 1043 : index = index - _num_sectors / 2 - 1;
375 1043 : ++counter;
376 : }
377 :
378 369 : counter = 0;
379 : // adding elements for other concentric circles
380 369 : index = Utility::pow<2>(standard / 2 + 1);
381 369 : limit = Utility::pow<2>(standard / 2 + 1) +
382 369 : (_num_sectors + 1) * (total_concentric_circles.size() - 1);
383 :
384 369 : int num_nodes_boundary = Utility::pow<2>(standard / 2 + 1) + standard + 1;
385 :
386 7147 : while (index < limit)
387 : {
388 6778 : Elem * elem = mesh->add_elem(std::make_unique<Quad4>());
389 6778 : elem->set_node(0, nodes[index]);
390 6778 : elem->set_node(1, nodes[index + standard + 1]);
391 6778 : elem->set_node(2, nodes[index + standard + 2]);
392 6778 : elem->set_node(3, nodes[index + 1]);
393 :
394 69940 : for (int i = 0; i < static_cast<int>(subdomainIDs.size() - 1); ++i)
395 63162 : if (index < limit - (standard + 1) * i && index >= limit - (standard + 1) * (i + 1))
396 6778 : elem->subdomain_id() = subdomainIDs[subdomainIDs.size() - 1 - i];
397 :
398 6778 : const int initial = Utility::pow<2>(standard / 2 + 1);
399 6778 : const int final = Utility::pow<2>(standard / 2 + 1) + standard - 1;
400 :
401 6778 : if ((index - initial) % (standard + 1) == 0)
402 1077 : boundary_info.add_side(elem, 0, 2);
403 6778 : if ((index - final) % (standard + 1) == 0)
404 1077 : boundary_info.add_side(elem, 2, 1);
405 6778 : if (!_has_outer_square)
406 3546 : if (index >= limit - (standard + 1))
407 694 : boundary_info.add_side(elem, 1, 3);
408 :
409 : // index increment is for adding nodes for a next element.
410 6778 : ++index;
411 :
412 : // increment by 2 indices may be necessary depending on where the index points to.
413 : // this varies based on the algorithms provided for the specific element and node placement.
414 6778 : if (index == (num_nodes_boundary + counter * (standard + 1)) - 1)
415 : {
416 1077 : ++index;
417 1077 : ++counter;
418 : }
419 : }
420 :
421 : // Enclosing square sections
422 : // ABCA
423 : // C B
424 : // B C
425 : // ACBA
426 :
427 : // adding elements for the enclosing square. (top left)
428 : int initial =
429 369 : Utility::pow<2>(standard / 2 + 1) + (standard + 1) * total_concentric_circles.size();
430 :
431 : int initial2 =
432 369 : Utility::pow<2>(standard / 2 + 1) + (standard + 1) * total_concentric_circles.size();
433 :
434 369 : if (_has_outer_square)
435 : {
436 104 : if (_rings.back() != 0) // this must be condition up front.
437 : {
438 104 : index = Utility::pow<2>(standard / 2 + 1) + (standard + 1) * total_concentric_circles.size();
439 104 : limit = Utility::pow<2>(standard / 2 + 1) + (standard + 1) * total_concentric_circles.size() +
440 104 : (_rings.back() + 1) * (standard / 2) + Utility::pow<2>(_rings.back() + 2) - 3 -
441 104 : _rings.back() * (_rings.back() + 2) - (_rings.back() + 1);
442 972 : while (index <= limit)
443 : {
444 : // outer square sector C
445 764 : Elem * elem = mesh->add_elem(std::make_unique<Quad4>());
446 764 : elem->set_node(0, nodes[index]);
447 764 : elem->set_node(1, nodes[index + 1]);
448 764 : elem->set_node(2, nodes[index + 1 + _rings.back() + 1]);
449 764 : elem->set_node(3, nodes[index + 1 + _rings.back()]);
450 764 : elem->subdomain_id() = subdomainIDs.back() + 1;
451 :
452 764 : if (index < (initial2 + static_cast<int>(_rings.back())))
453 206 : boundary_info.add_side(elem, 0, 1);
454 :
455 764 : if (index == initial)
456 356 : boundary_info.add_side(elem, 3, 4);
457 :
458 764 : ++index;
459 :
460 : // As mentioned before, increment by 2 indices may be required depending on where the index
461 : // points to.
462 764 : if ((index - initial) % static_cast<int>(_rings.back()) == 0)
463 : {
464 356 : ++index;
465 356 : initial = initial + (static_cast<int>(_rings.back()) + 1);
466 : }
467 : }
468 :
469 : // adding elements for the enclosing square. (top right)
470 104 : initial = Utility::pow<2>(standard / 2 + 1) +
471 312 : (standard + 1) * total_concentric_circles.size() +
472 104 : (_rings.back() + 1) * (standard / 2 + 1);
473 :
474 104 : limit = Utility::pow<2>(standard / 2 + 1) + (standard + 1) * total_concentric_circles.size() +
475 104 : (_rings.back() + 1) * (standard / 2) + Utility::pow<2>(_rings.back() + 2) - 3 -
476 104 : (_rings.back() + 2);
477 :
478 1016 : while (index <= limit)
479 : {
480 : // outer square sector A
481 808 : Elem * elem = mesh->add_elem(std::make_unique<Quad4>());
482 808 : elem->set_node(3, nodes[index]);
483 808 : elem->set_node(2, nodes[index + _rings.back() + 2]);
484 808 : elem->set_node(1, nodes[index + _rings.back() + 3]);
485 808 : elem->set_node(0, nodes[index + 1]);
486 808 : elem->subdomain_id() = subdomainIDs.back() + 1;
487 :
488 808 : if (index >= static_cast<int>(limit - (_rings.back() + 1)))
489 310 : boundary_info.add_side(elem, 1, 3);
490 :
491 808 : if ((index - initial) % static_cast<int>(_rings.back() + 2) == 0)
492 206 : boundary_info.add_side(elem, 2, 4);
493 :
494 808 : ++index;
495 :
496 808 : if ((index - initial) % static_cast<int>(_rings.back() + 1) == 0)
497 : {
498 206 : ++index;
499 206 : initial = initial + (static_cast<int>(_rings.back()) + 2);
500 : }
501 : }
502 :
503 : // adding elements for the enclosing square. (one center quad)
504 104 : int index1 = Utility::pow<2>(standard / 2 + 1) +
505 104 : (standard + 1) * (total_concentric_circles.size() - 1) + standard / 2;
506 :
507 104 : int index2 = Utility::pow<2>(standard / 2 + 1) +
508 312 : (standard + 1) * total_concentric_circles.size() +
509 104 : (_rings.back() + 1) * (standard / 2) + Utility::pow<2>(_rings.back() + 2) - 3 -
510 104 : _rings.back() * (_rings.back() + 2) - (_rings.back() + 1);
511 :
512 : // pointy tips of the A sectors, touching the inner circle
513 104 : Elem * elem = mesh->add_elem(std::make_unique<Quad4>());
514 104 : elem->set_node(3, nodes[index1]);
515 104 : elem->set_node(2, nodes[index2]);
516 104 : elem->set_node(1, nodes[index2 + _rings.back() + 1]);
517 104 : elem->set_node(0, nodes[index2 + _rings.back() + 2]);
518 104 : elem->subdomain_id() = subdomainIDs.back() + 1;
519 :
520 : // adding elements for the left mid part.
521 104 : index = Utility::pow<2>(standard / 2 + 1) + standard / 2 +
522 104 : (standard + 1) * (total_concentric_circles.size() - 1);
523 104 : limit = index + standard / 2 - 1;
524 :
525 356 : while (index <= limit)
526 : {
527 : // outer square elements in sector C touching the inner circle
528 252 : Elem * elem = mesh->add_elem(std::make_unique<Quad4>());
529 252 : elem->set_node(3, nodes[index]);
530 252 : elem->set_node(2, nodes[index + 1]);
531 252 : elem->set_node(1, nodes[index2 - _rings.back() - 1]);
532 252 : elem->set_node(0, nodes[index2]);
533 252 : elem->subdomain_id() = subdomainIDs.back() + 1;
534 :
535 252 : if (index == limit)
536 104 : boundary_info.add_side(elem, 1, 1);
537 :
538 252 : ++index;
539 :
540 : // two different indices are used to add nodes for an element.
541 252 : index2 = index2 - _rings.back() - 1;
542 : }
543 :
544 : // adding elements for the right mid part.
545 104 : index1 = Utility::pow<2>(standard / 2 + 1) + standard / 2 +
546 104 : (standard + 1) * (total_concentric_circles.size() - 1);
547 104 : index2 = Utility::pow<2>(standard / 2 + 1) +
548 312 : (standard + 1) * total_concentric_circles.size() +
549 104 : (_rings.back() + 1) * (standard / 2) + Utility::pow<2>(_rings.back() + 2) - 2 +
550 104 : (_rings.back() + 1);
551 : int index3 =
552 104 : Utility::pow<2>(standard / 2 + 1) + (standard + 1) * total_concentric_circles.size() +
553 104 : (_rings.back() + 1) * (standard / 2) - 1 + (_rings.back() + 1) + (_rings.back() + 2);
554 :
555 : // elements clockwise from the A sector tips
556 104 : elem = mesh->add_elem(std::make_unique<Quad4>());
557 104 : elem->set_node(0, nodes[index1]);
558 104 : elem->set_node(1, nodes[index1 - 1]);
559 104 : elem->set_node(2, nodes[index2]);
560 104 : elem->set_node(3, nodes[index3]);
561 104 : elem->subdomain_id() = subdomainIDs.back() + 1;
562 :
563 104 : if (standard == 2)
564 10 : boundary_info.add_side(elem, 1, 2);
565 :
566 : // adding elements for the right mid bottom part.
567 :
568 104 : index = Utility::pow<2>(standard / 2 + 1) + standard / 2 +
569 104 : (standard + 1) * (total_concentric_circles.size() - 1) - 2;
570 104 : index1 = Utility::pow<2>(standard / 2 + 1) +
571 312 : (standard + 1) * total_concentric_circles.size() +
572 104 : (_rings.back() + 1) * (standard / 2) + Utility::pow<2>(_rings.back() + 2) - 2 +
573 104 : (_rings.back() + 1) * 2;
574 :
575 104 : limit = Utility::pow<2>(standard / 2 + 1) + standard / 2 +
576 104 : (standard + 1) * (total_concentric_circles.size() - 1) - standard / 2;
577 :
578 104 : if (standard != 2)
579 : {
580 242 : while (index >= limit)
581 : {
582 : // outer square elements in sector B touching the inner circle
583 148 : Elem * elem = mesh->add_elem(std::make_unique<Quad4>());
584 148 : elem->set_node(0, nodes[index]);
585 148 : elem->set_node(1, nodes[index1]);
586 148 : elem->set_node(2, nodes[index1 - (_rings.back() + 1)]);
587 148 : elem->set_node(3, nodes[index + 1]);
588 148 : elem->subdomain_id() = subdomainIDs.back() + 1;
589 :
590 148 : if (index == limit)
591 94 : boundary_info.add_side(elem, 0, 2);
592 148 : --index;
593 148 : index1 = index1 + (_rings.back() + 1);
594 : }
595 : }
596 :
597 : // adding elements for the right low part.
598 104 : index = Utility::pow<2>(standard / 2 + 1) + (standard + 1) * total_concentric_circles.size() +
599 104 : (_rings.back() + 1) * (standard / 2) + Utility::pow<2>(_rings.back() + 2) - 2;
600 :
601 104 : index1 = index - (_rings.back() + 2);
602 : // dummy condition for elem definition
603 104 : if (standard >= 2)
604 : {
605 : // single elements between A and B on the outside of the square
606 104 : Elem * elem = mesh->add_elem(std::make_unique<Quad4>());
607 104 : elem->set_node(3, nodes[index]);
608 104 : elem->set_node(2, nodes[index + 1]);
609 104 : elem->set_node(1, nodes[index + 2]);
610 104 : elem->set_node(0, nodes[index1]);
611 104 : elem->subdomain_id() = subdomainIDs.back() + 1;
612 :
613 104 : boundary_info.add_side(elem, 2, 3);
614 :
615 104 : if (standard == 2)
616 10 : boundary_info.add_side(elem, 1, 2);
617 : }
618 :
619 104 : index = Utility::pow<2>(standard / 2 + 1) + (standard + 1) * total_concentric_circles.size() +
620 104 : (_rings.back() + 1) * (standard / 2) + Utility::pow<2>(_rings.back() + 2) - 2 -
621 104 : (_rings.back() + 2);
622 :
623 104 : limit = Utility::pow<2>(standard / 2 + 1) + (standard + 1) * total_concentric_circles.size() +
624 104 : (_rings.back() + 1) * (standard / 2) + (_rings.back() + 2) * 2 - 2;
625 :
626 104 : int k = 1;
627 206 : while (index > limit)
628 : {
629 102 : Elem * elem = mesh->add_elem(std::make_unique<Quad4>());
630 102 : elem->set_node(3, nodes[index]);
631 102 : elem->set_node(2, nodes[index + (_rings.back() + 2) * k + k + 1]);
632 102 : elem->set_node(1, nodes[index + (_rings.back() + 2) * k + k + 2]);
633 102 : elem->set_node(0, nodes[index - _rings.back() - 2]);
634 102 : elem->subdomain_id() = subdomainIDs.back() + 1;
635 102 : index = index - (_rings.back() + 2);
636 102 : ++k;
637 :
638 102 : if (standard == 2)
639 0 : boundary_info.add_side(elem, 1, 2);
640 : }
641 :
642 104 : index = Utility::pow<2>(standard / 2 + 1) + (standard + 1) * total_concentric_circles.size() +
643 104 : (_rings.back() + 1) * (standard / 2) + Utility::pow<2>(_rings.back() + 2) - 1;
644 104 : initial = Utility::pow<2>(standard / 2 + 1) +
645 312 : (standard + 1) * total_concentric_circles.size() +
646 104 : (_rings.back() + 1) * (standard / 2) + Utility::pow<2>(_rings.back() + 2) - 1;
647 104 : limit = Utility::pow<2>(standard / 2 + 1) + (standard + 1) * total_concentric_circles.size() +
648 104 : (_rings.back() + 1) * (standard / 2) * 2 + Utility::pow<2>(_rings.back() + 2) - 2 -
649 104 : _rings.back() - 1;
650 :
651 104 : initial2 = Utility::pow<2>(standard / 2 + 1) +
652 312 : (standard + 1) * total_concentric_circles.size() +
653 104 : (_rings.back() + 1) * (standard / 2) * 2 + Utility::pow<2>(_rings.back() + 2) - 2 -
654 104 : (_rings.back() + 1) * 2;
655 :
656 104 : if (standard > 2)
657 : {
658 540 : while (index < limit)
659 : {
660 352 : Elem * elem = mesh->add_elem(std::make_unique<Quad4>());
661 352 : elem->set_node(0, nodes[index]);
662 352 : elem->set_node(1, nodes[index + 1]);
663 352 : elem->set_node(2, nodes[index + 1 + _rings.back() + 1]);
664 352 : elem->set_node(3, nodes[index + 1 + _rings.back()]);
665 352 : elem->subdomain_id() = subdomainIDs.back() + 1;
666 :
667 352 : if (index > initial2)
668 196 : boundary_info.add_side(elem, 2, 2);
669 :
670 352 : if ((index - initial) == 0)
671 148 : boundary_info.add_side(elem, 3, 3);
672 :
673 352 : ++index;
674 :
675 352 : if ((index - initial) % static_cast<int>(_rings.back()) == 0)
676 : {
677 148 : ++index;
678 148 : initial = initial + (static_cast<int>(_rings.back()) + 1);
679 : }
680 : }
681 : }
682 : }
683 : }
684 :
685 : // This is to set boundary names.
686 369 : boundary_info.sideset_name(1) = "left";
687 369 : boundary_info.sideset_name(2) = "bottom";
688 :
689 369 : if (!_has_outer_square)
690 265 : boundary_info.sideset_name(3) = "outer";
691 : else
692 : {
693 104 : boundary_info.sideset_name(3) = "right";
694 104 : boundary_info.sideset_name(4) = "top";
695 : }
696 :
697 369 : if (_portion == "top_left")
698 : {
699 10 : MeshTools::Modification::rotate(*mesh, 90, 0, 0);
700 10 : boundary_info.sideset_name(1) = "bottom";
701 10 : boundary_info.sideset_name(2) = "right";
702 :
703 10 : if (!_has_outer_square)
704 10 : boundary_info.sideset_name(3) = "outer";
705 : else
706 : {
707 0 : boundary_info.sideset_name(3) = "top";
708 0 : boundary_info.sideset_name(4) = "left";
709 : }
710 : }
711 359 : else if (_portion == "bottom_left")
712 : {
713 10 : MeshTools::Modification::rotate(*mesh, 180, 0, 0);
714 10 : boundary_info.sideset_name(1) = "right";
715 10 : boundary_info.sideset_name(2) = "top";
716 :
717 10 : if (!_has_outer_square)
718 10 : boundary_info.sideset_name(3) = "outer";
719 : else
720 : {
721 0 : boundary_info.sideset_name(3) = "left";
722 0 : boundary_info.sideset_name(4) = "bottom";
723 : }
724 : }
725 349 : else if (_portion == "bottom_right")
726 : {
727 10 : MeshTools::Modification::rotate(*mesh, 270, 0, 0);
728 10 : boundary_info.sideset_name(1) = "top";
729 10 : boundary_info.sideset_name(2) = "left";
730 :
731 10 : if (!_has_outer_square)
732 10 : boundary_info.sideset_name(3) = "outer";
733 : else
734 : {
735 0 : boundary_info.sideset_name(3) = "bottom";
736 0 : boundary_info.sideset_name(4) = "right";
737 : }
738 : }
739 :
740 339 : else if (_portion == "top_half")
741 : {
742 0 : ReplicatedMesh other_mesh(*mesh);
743 : // This is to rotate the mesh and also to reset boundary IDs.
744 0 : MeshTools::Modification::rotate(other_mesh, 90, 0, 0);
745 0 : if (_has_outer_square)
746 : {
747 0 : MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
748 0 : MeshTools::Modification::change_boundary_id(other_mesh, 2, 6);
749 0 : MeshTools::Modification::change_boundary_id(other_mesh, 3, 7);
750 0 : MeshTools::Modification::change_boundary_id(other_mesh, 4, 1);
751 0 : MeshTools::Modification::change_boundary_id(other_mesh, 5, 2);
752 0 : MeshTools::Modification::change_boundary_id(other_mesh, 6, 3);
753 0 : MeshTools::Modification::change_boundary_id(other_mesh, 7, 4);
754 0 : mesh->prepare_for_use();
755 0 : other_mesh.prepare_for_use();
756 0 : mesh->stitch_meshes(other_mesh, 1, 3, TOLERANCE, true, /*verbose=*/false);
757 0 : mesh->get_boundary_info().sideset_name(1) = "left";
758 0 : mesh->get_boundary_info().sideset_name(2) = "bottom";
759 0 : mesh->get_boundary_info().sideset_name(3) = "right";
760 0 : mesh->get_boundary_info().sideset_name(4) = "top";
761 : }
762 : else
763 : {
764 0 : MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
765 0 : MeshTools::Modification::change_boundary_id(other_mesh, 2, 1);
766 0 : MeshTools::Modification::change_boundary_id(other_mesh, 5, 2);
767 0 : mesh->prepare_for_use();
768 0 : other_mesh.prepare_for_use();
769 0 : mesh->stitch_meshes(other_mesh, 1, 1, TOLERANCE, true, /*verbose=*/false);
770 :
771 0 : MeshTools::Modification::change_boundary_id(*mesh, 2, 1);
772 0 : MeshTools::Modification::change_boundary_id(*mesh, 3, 2);
773 0 : mesh->get_boundary_info().sideset_name(1) = "bottom";
774 0 : mesh->get_boundary_info().sideset_name(2) = "outer";
775 : }
776 0 : other_mesh.clear();
777 0 : }
778 :
779 339 : else if (_portion == "right_half")
780 : {
781 0 : ReplicatedMesh other_mesh(*mesh);
782 : // This is to rotate the mesh and also to reset boundary IDs.
783 0 : MeshTools::Modification::rotate(other_mesh, 270, 0, 0);
784 0 : if (_has_outer_square)
785 : {
786 0 : MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
787 0 : MeshTools::Modification::change_boundary_id(other_mesh, 2, 6);
788 0 : MeshTools::Modification::change_boundary_id(other_mesh, 3, 7);
789 0 : MeshTools::Modification::change_boundary_id(other_mesh, 4, 3);
790 0 : MeshTools::Modification::change_boundary_id(other_mesh, 5, 4);
791 0 : MeshTools::Modification::change_boundary_id(other_mesh, 6, 1);
792 0 : MeshTools::Modification::change_boundary_id(other_mesh, 7, 2);
793 0 : mesh->prepare_for_use();
794 0 : other_mesh.prepare_for_use();
795 0 : mesh->stitch_meshes(other_mesh, 2, 4, TOLERANCE, true, /*verbose=*/false);
796 0 : mesh->get_boundary_info().sideset_name(1) = "left";
797 0 : mesh->get_boundary_info().sideset_name(2) = "bottom";
798 0 : mesh->get_boundary_info().sideset_name(3) = "right";
799 0 : mesh->get_boundary_info().sideset_name(4) = "top";
800 : }
801 : else
802 : {
803 0 : MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
804 0 : MeshTools::Modification::change_boundary_id(other_mesh, 2, 1);
805 0 : MeshTools::Modification::change_boundary_id(other_mesh, 5, 2);
806 0 : mesh->prepare_for_use();
807 0 : other_mesh.prepare_for_use();
808 0 : mesh->stitch_meshes(other_mesh, 2, 2, TOLERANCE, true, /*verbose=*/false);
809 :
810 0 : MeshTools::Modification::change_boundary_id(*mesh, 3, 2);
811 0 : mesh->get_boundary_info().sideset_name(1) = "left";
812 0 : mesh->get_boundary_info().sideset_name(2) = "outer";
813 : }
814 0 : other_mesh.clear();
815 0 : }
816 339 : else if (_portion == "left_half")
817 : {
818 0 : ReplicatedMesh other_mesh(*mesh);
819 :
820 : // This is to rotate the mesh and to reset boundary IDs.
821 0 : MeshTools::Modification::rotate(other_mesh, 90, 0, 0);
822 0 : MeshTools::Modification::rotate(*mesh, 180, 0, 0);
823 0 : if (_has_outer_square)
824 : {
825 : // The other mesh is created by rotating the original mesh about 90 degrees.
826 0 : MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
827 0 : MeshTools::Modification::change_boundary_id(other_mesh, 2, 6);
828 0 : MeshTools::Modification::change_boundary_id(other_mesh, 3, 7);
829 0 : MeshTools::Modification::change_boundary_id(other_mesh, 4, 1);
830 0 : MeshTools::Modification::change_boundary_id(other_mesh, 5, 2);
831 0 : MeshTools::Modification::change_boundary_id(other_mesh, 6, 3);
832 0 : MeshTools::Modification::change_boundary_id(other_mesh, 7, 4);
833 : // The original mesh is then rotated about 180 degrees.
834 0 : MeshTools::Modification::change_boundary_id(*mesh, 1, 5);
835 0 : MeshTools::Modification::change_boundary_id(*mesh, 2, 6);
836 0 : MeshTools::Modification::change_boundary_id(*mesh, 3, 7);
837 0 : MeshTools::Modification::change_boundary_id(*mesh, 4, 2);
838 0 : MeshTools::Modification::change_boundary_id(*mesh, 5, 3);
839 0 : MeshTools::Modification::change_boundary_id(*mesh, 6, 4);
840 0 : MeshTools::Modification::change_boundary_id(*mesh, 7, 1);
841 0 : mesh->prepare_for_use();
842 0 : other_mesh.prepare_for_use();
843 0 : mesh->stitch_meshes(other_mesh, 4, 2, TOLERANCE, true, /*verbose=*/false);
844 0 : mesh->get_boundary_info().sideset_name(1) = "left";
845 0 : mesh->get_boundary_info().sideset_name(2) = "bottom";
846 0 : mesh->get_boundary_info().sideset_name(3) = "right";
847 0 : mesh->get_boundary_info().sideset_name(4) = "top";
848 : }
849 : else
850 : {
851 0 : MeshTools::Modification::change_boundary_id(*mesh, 1, 5);
852 0 : MeshTools::Modification::change_boundary_id(*mesh, 2, 1);
853 0 : MeshTools::Modification::change_boundary_id(*mesh, 5, 2);
854 0 : mesh->prepare_for_use();
855 0 : other_mesh.prepare_for_use();
856 0 : mesh->stitch_meshes(other_mesh, 1, 1, TOLERANCE, true, /*verbose=*/false);
857 :
858 0 : MeshTools::Modification::change_boundary_id(*mesh, 2, 1);
859 0 : MeshTools::Modification::change_boundary_id(*mesh, 3, 2);
860 0 : mesh->get_boundary_info().sideset_name(1) = "right";
861 0 : mesh->get_boundary_info().sideset_name(2) = "outer";
862 : }
863 0 : other_mesh.clear();
864 0 : }
865 339 : else if (_portion == "bottom_half")
866 : {
867 0 : ReplicatedMesh other_mesh(*mesh);
868 : // This is to rotate the mesh and also to reset boundary IDs.
869 0 : MeshTools::Modification::rotate(other_mesh, 180, 0, 0);
870 0 : MeshTools::Modification::rotate(*mesh, 270, 0, 0);
871 0 : if (_has_outer_square)
872 : {
873 : // The other mesh is created by rotating the original mesh about 180 degrees.
874 0 : MeshTools::Modification::change_boundary_id(other_mesh, 1, 5);
875 0 : MeshTools::Modification::change_boundary_id(other_mesh, 2, 6);
876 0 : MeshTools::Modification::change_boundary_id(other_mesh, 3, 7);
877 0 : MeshTools::Modification::change_boundary_id(other_mesh, 4, 2);
878 0 : MeshTools::Modification::change_boundary_id(other_mesh, 5, 3);
879 0 : MeshTools::Modification::change_boundary_id(other_mesh, 6, 4);
880 0 : MeshTools::Modification::change_boundary_id(other_mesh, 7, 1);
881 : // The original mesh is rotated about 270 degrees.
882 0 : MeshTools::Modification::change_boundary_id(*mesh, 1, 5);
883 0 : MeshTools::Modification::change_boundary_id(*mesh, 2, 6);
884 0 : MeshTools::Modification::change_boundary_id(*mesh, 3, 7);
885 0 : MeshTools::Modification::change_boundary_id(*mesh, 4, 3);
886 0 : MeshTools::Modification::change_boundary_id(*mesh, 5, 4);
887 0 : MeshTools::Modification::change_boundary_id(*mesh, 6, 1);
888 0 : MeshTools::Modification::change_boundary_id(*mesh, 7, 2);
889 0 : mesh->prepare_for_use();
890 0 : other_mesh.prepare_for_use();
891 0 : mesh->stitch_meshes(other_mesh, 1, 3, TOLERANCE, true, /*verbose=*/false);
892 0 : mesh->get_boundary_info().sideset_name(1) = "left";
893 0 : mesh->get_boundary_info().sideset_name(2) = "bottom";
894 0 : mesh->get_boundary_info().sideset_name(3) = "right";
895 0 : mesh->get_boundary_info().sideset_name(4) = "top";
896 : }
897 : else
898 : {
899 0 : MeshTools::Modification::change_boundary_id(*mesh, 1, 5);
900 0 : MeshTools::Modification::change_boundary_id(*mesh, 2, 1);
901 0 : MeshTools::Modification::change_boundary_id(*mesh, 5, 2);
902 0 : mesh->prepare_for_use();
903 0 : other_mesh.prepare_for_use();
904 0 : mesh->stitch_meshes(other_mesh, 1, 1, TOLERANCE, true, /*verbose=*/false);
905 :
906 0 : MeshTools::Modification::change_boundary_id(*mesh, 2, 1);
907 0 : MeshTools::Modification::change_boundary_id(*mesh, 3, 2);
908 0 : mesh->get_boundary_info().sideset_name(1) = "top";
909 0 : mesh->get_boundary_info().sideset_name(2) = "outer";
910 : }
911 0 : other_mesh.clear();
912 0 : }
913 339 : else if (_portion == "full")
914 : {
915 249 : ReplicatedMesh portion_two(*mesh);
916 :
917 : // This is to rotate the mesh and also to reset boundary IDs.
918 249 : MeshTools::Modification::rotate(portion_two, 90, 0, 0);
919 :
920 249 : if (_has_outer_square)
921 : {
922 : // Portion 2: 2nd quadrant
923 64 : MeshTools::Modification::change_boundary_id(portion_two, 1, 5);
924 64 : MeshTools::Modification::change_boundary_id(portion_two, 2, 6);
925 64 : MeshTools::Modification::change_boundary_id(portion_two, 3, 7);
926 64 : MeshTools::Modification::change_boundary_id(portion_two, 4, 1);
927 64 : MeshTools::Modification::change_boundary_id(portion_two, 5, 2);
928 64 : MeshTools::Modification::change_boundary_id(portion_two, 6, 3);
929 64 : MeshTools::Modification::change_boundary_id(portion_two, 7, 4);
930 64 : mesh->prepare_for_use();
931 64 : portion_two.prepare_for_use();
932 : // 'top_half'
933 64 : mesh->stitch_meshes(portion_two, 1, 3, TOLERANCE, true, /*verbose=*/false);
934 :
935 : // 'bottom_half'
936 64 : ReplicatedMesh portion_bottom(*mesh);
937 64 : MeshTools::Modification::rotate(portion_bottom, 180, 0, 0);
938 64 : MeshTools::Modification::change_boundary_id(portion_bottom, 1, 5);
939 64 : MeshTools::Modification::change_boundary_id(portion_bottom, 2, 6);
940 64 : MeshTools::Modification::change_boundary_id(portion_bottom, 3, 7);
941 64 : MeshTools::Modification::change_boundary_id(portion_bottom, 4, 2);
942 64 : MeshTools::Modification::change_boundary_id(portion_bottom, 5, 3);
943 64 : MeshTools::Modification::change_boundary_id(portion_bottom, 6, 4);
944 64 : MeshTools::Modification::change_boundary_id(portion_bottom, 7, 1);
945 64 : mesh->prepare_for_use();
946 64 : portion_bottom.prepare_for_use();
947 : // 'full'
948 64 : mesh->stitch_meshes(portion_bottom, 2, 4, TOLERANCE, true, /*verbose=*/false);
949 :
950 64 : mesh->get_boundary_info().sideset_name(1) = "left";
951 64 : mesh->get_boundary_info().sideset_name(2) = "bottom";
952 64 : mesh->get_boundary_info().sideset_name(3) = "right";
953 64 : mesh->get_boundary_info().sideset_name(4) = "top";
954 64 : portion_bottom.clear();
955 64 : }
956 : else
957 : {
958 185 : MeshTools::Modification::change_boundary_id(portion_two, 1, 5);
959 185 : MeshTools::Modification::change_boundary_id(portion_two, 2, 1);
960 185 : MeshTools::Modification::change_boundary_id(portion_two, 5, 2);
961 : // 'top half'
962 185 : mesh->prepare_for_use();
963 185 : portion_two.prepare_for_use();
964 185 : mesh->stitch_meshes(portion_two, 1, 1, TOLERANCE, true, /*verbose=*/false);
965 : // 'bottom half'
966 185 : ReplicatedMesh portion_bottom(*mesh);
967 185 : MeshTools::Modification::rotate(portion_bottom, 180, 0, 0);
968 : // 'full'
969 185 : mesh->prepare_for_use();
970 185 : portion_bottom.prepare_for_use();
971 185 : mesh->stitch_meshes(portion_bottom, 2, 2, TOLERANCE, true, /*verbose=*/false);
972 185 : MeshTools::Modification::change_boundary_id(*mesh, 3, 1);
973 185 : mesh->get_boundary_info().sideset_name(1) = "outer";
974 185 : portion_bottom.clear();
975 185 : }
976 249 : portion_two.clear();
977 249 : }
978 :
979 738 : if (_portion != "top_half" && _portion != "right_half" && _portion != "left_half" &&
980 738 : _portion != "bottom_half" && _portion != "full")
981 120 : mesh->prepare_for_use();
982 :
983 : // Laplace smoothing
984 369 : libMesh::LaplaceMeshSmoother lms(*mesh, _smoothing_max_it);
985 369 : lms.smooth();
986 :
987 : // Add subdomain name
988 1107 : if (isParamValid("subdomain_name"))
989 0 : mesh->subdomain_name(0) = getParam<SubdomainName>("subdomain_name");
990 :
991 369 : mesh->prepare_for_use();
992 738 : return dynamic_pointer_cast<MeshBase>(mesh);
993 369 : }
|
