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