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