www.mooseframework.org
SpiralAnnularMeshGenerator.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 
13 #include "libmesh/replicated_mesh.h"
14 #include "libmesh/face_quad4.h"
15 #include "libmesh/face_tri3.h"
16 
18 
19 template <>
22 {
24 
25  params.addRequiredRangeCheckedParam<Real>(
26  "inner_radius", "inner_radius>0.", "The size of the inner circle.");
27  params.addRequiredRangeCheckedParam<Real>("outer_radius",
28  "outer_radius>0.",
29  "The size of the outer circle."
30  " Logically, it has to be greater than inner_radius");
31  params.addRequiredRangeCheckedParam<unsigned int>(
32  "nodes_per_ring", "nodes_per_ring>5", "Number of nodes on each ring.");
33  params.addParam<bool>(
34  "use_tri6", false, "Generate mesh of TRI6 elements instead of TRI3 elements.");
35  params.addRequiredRangeCheckedParam<unsigned int>(
36  "num_rings", "num_rings>1", "The number of rings.");
37  params.addParam<boundary_id_type>(
38  "cylinder_bid", 1, "The boundary id to use for the cylinder (inner circle)");
39  params.addParam<boundary_id_type>(
40  "exterior_bid", 2, "The boundary id to use for the exterior (outer circle)");
41  params.addParam<Real>("initial_delta_r",
42  "Width of the initial layer of elements around the cylinder."
43  "This number should be approximately"
44  " 2 * pi * inner_radius / nodes_per_ring to ensure that the"
45  " initial layer of elements is almost equilateral");
46  params.addClassDescription("Creates an annual mesh based on TRI3 elements"
47  " (it can also be TRI6 elements) on several rings.");
48 
49  return params;
50 }
51 
53  : MeshGenerator(parameters),
54  _inner_radius(getParam<Real>("inner_radius")),
55  _outer_radius(getParam<Real>("outer_radius")),
56  _radial_bias(1.0),
57  _nodes_per_ring(getParam<unsigned int>("nodes_per_ring")),
58  _use_tri6(getParam<bool>("use_tri6")),
59  _num_rings(getParam<unsigned int>("num_rings")),
60  _cylinder_bid(getParam<boundary_id_type>("cylinder_bid")),
61  _exterior_bid(getParam<boundary_id_type>("exterior_bid")),
62  _initial_delta_r(2 * libMesh::pi * _inner_radius / _nodes_per_ring)
63 {
64  // catch likely user errors
66  mooseError("SpiralAnnularMesh: outer_radius must be greater than inner_radius");
67 }
68 
69 std::unique_ptr<MeshBase>
71 {
72  std::unique_ptr<ReplicatedMesh> mesh = libmesh_make_unique<ReplicatedMesh>(comm(), 2);
73 
74  {
75  // Compute the radial bias given:
76  // .) the inner radius
77  // .) the outer radius
78  // .) the initial_delta_r
79  // .) the desired number of intervals
80  // Note: the exponent n used in the formula is one less than the
81  // number of rings the user requests.
82  Real alpha = 1.1;
83  int n = _num_rings - 1;
84 
85  // lambda used to compute the residual and Jacobian for the Newton iterations.
86  // We capture parameters which don't need to change from the current scope at
87  // the time this lambda is declared. The values are not updated later, so we
88  // can't use this for e.g. f, df, and alpha.
89  auto newton = [this, n](Real & f, Real & df, const Real & alpha) {
90  f = (1. - std::pow(alpha, n + 1)) / (1. - alpha) -
92  df = (-(n + 1) * (1 - alpha) * std::pow(alpha, n) + (1. - std::pow(alpha, n + 1))) /
93  (1. - alpha) / (1. - alpha);
94  };
95 
96  Real f, df;
97  int num_iter = 1;
98  newton(f, df, alpha);
99 
100  while (std::abs(f) > 1.e-9 && num_iter <= 25)
101  {
102  // Compute and apply update.
103  Real dx = -f / df;
104  alpha += dx;
105  newton(f, df, alpha);
106  num_iter++;
107  }
108 
109  // In case the Newton iteration fails to converge.
110  if (num_iter > 25)
111  mooseError("Newton iteration failed to converge (more than 25 iterations).");
112 
113  // Set radial basis to the value of alpha that we computed with Newton.
114  _radial_bias = alpha;
115  }
116 
117  // The number of rings specified by the user does not include the ring at
118  // the surface of the cylinder itself, so we increment it by one now.
119  _num_rings += 1;
120 
121  // Data structure that holds pointers to the Nodes of each ring.
122  std::vector<std::vector<Node *>> ring_nodes(_num_rings);
123 
124  // Initialize radius and delta_r variables.
125  Real radius = _inner_radius;
126  Real delta_r = _initial_delta_r;
127 
128  // Node id counter.
129  unsigned int current_node_id = 0;
130 
131  for (std::size_t r = 0; r < _num_rings; ++r)
132  {
133  ring_nodes[r].resize(_nodes_per_ring);
134 
135  // Add nodes starting from either theta=0 or theta=pi/nodes_per_ring
136  Real theta = r % 2 == 0 ? 0 : (libMesh::pi / _nodes_per_ring);
137  for (std::size_t n = 0; n < _nodes_per_ring; ++n)
138  {
139  ring_nodes[r][n] = mesh->add_point(Point(radius * std::cos(theta), radius * std::sin(theta)),
140  current_node_id++);
141  // Update angle
142  theta += 2 * libMesh::pi / _nodes_per_ring;
143  }
144 
145  // Go to next ring
146  radius += delta_r;
147  delta_r *= _radial_bias;
148  }
149 
150  // Add elements
151  for (std::size_t r = 0; r < _num_rings - 1; ++r)
152  {
153  // even -> odd ring
154  if (r % 2 == 0)
155  {
156  // Inner ring (n, n*, n+1)
157  // Starred indices refer to nodes on the "outer" ring of this pair.
158  for (std::size_t n = 0; n < _nodes_per_ring; ++n)
159  {
160  // Wrap around
161  unsigned int np1 = (n == _nodes_per_ring - 1) ? 0 : n + 1;
162  Elem * elem = mesh->add_elem(new Tri3);
163  elem->set_node(0) = ring_nodes[r][n];
164  elem->set_node(1) = ring_nodes[r + 1][n];
165  elem->set_node(2) = ring_nodes[r][np1];
166 
167  // Add interior faces to 'cylinder' sideset if we are on ring 0.
168  if (r == 0)
169  mesh->boundary_info->add_side(elem->id(), /*side=*/2, _cylinder_bid);
170  }
171 
172  // Outer ring (n*, n+1*, n+1)
173  for (std::size_t n = 0; n < _nodes_per_ring; ++n)
174  {
175  // Wrap around
176  unsigned int np1 = (n == _nodes_per_ring - 1) ? 0 : n + 1;
177  Elem * elem = mesh->add_elem(new Tri3);
178  elem->set_node(0) = ring_nodes[r + 1][n];
179  elem->set_node(1) = ring_nodes[r + 1][np1];
180  elem->set_node(2) = ring_nodes[r][np1];
181 
182  // Add exterior faces to 'exterior' sideset if we're on the last ring.
183  // Note: this code appears in two places since we could end on either an even or odd ring.
184  if (r == _num_rings - 2)
185  mesh->boundary_info->add_side(elem->id(), /*side=*/0, _exterior_bid);
186  }
187  }
188  else
189  {
190  // odd -> even ring
191  // Inner ring (n, n+1*, n+1)
192  for (std::size_t n = 0; n < _nodes_per_ring; ++n)
193  {
194  // Wrap around
195  unsigned int np1 = (n == _nodes_per_ring - 1) ? 0 : n + 1;
196  Elem * elem = mesh->add_elem(new Tri3);
197  elem->set_node(0) = ring_nodes[r][n];
198  elem->set_node(1) = ring_nodes[r + 1][np1];
199  elem->set_node(2) = ring_nodes[r][np1];
200  }
201 
202  // Outer ring (n*, n+1*, n)
203  for (std::size_t n = 0; n < _nodes_per_ring; ++n)
204  {
205  // Wrap around
206  unsigned int np1 = (n == _nodes_per_ring - 1) ? 0 : n + 1;
207  Elem * elem = mesh->add_elem(new Tri3);
208  elem->set_node(0) = ring_nodes[r + 1][n];
209  elem->set_node(1) = ring_nodes[r + 1][np1];
210  elem->set_node(2) = ring_nodes[r][n];
211 
212  // Add exterior faces to 'exterior' sideset if we're on the last ring.
213  if (r == _num_rings - 2)
214  mesh->boundary_info->add_side(elem->id(), /*side=*/0, _exterior_bid);
215  }
216  }
217  }
218 
219  // Sanity check: make sure all elements have positive area. Note: we
220  // can't use elem->volume() for this, as that always returns a
221  // positive area regardless of the node ordering.
222  // We compute (p1-p0) \cross (p2-p0) and check that the z-component is positive.
223  for (const auto & elem : mesh->element_ptr_range())
224  {
225  Point cp = (elem->point(1) - elem->point(0)).cross(elem->point(2) - elem->point(0));
226  if (cp(2) < 0.)
227  mooseError("Invalid elem found with negative area");
228  }
229 
230  // Create sideset names.
231  mesh->boundary_info->sideset_name(_cylinder_bid) = "cylinder";
232  mesh->boundary_info->sideset_name(_exterior_bid) = "exterior";
233 
234  // Find neighbors, etc.
235  mesh->prepare_for_use();
236 
237  if (_use_tri6)
238  {
239  mesh->all_second_order(/*full_ordered=*/true);
240  std::vector<unsigned int> nos;
241 
242  // Loop over the elements, moving mid-edge nodes onto the
243  // nearest radius as applicable. For each element, exactly one
244  // edge should lie on the same radius, so we move only that
245  // mid-edge node.
246  for (const auto & elem : mesh->element_ptr_range())
247  {
248  // Make sure we are dealing only with triangles
249  libmesh_assert(elem->n_vertices() == 3);
250 
251  // Compute vertex radii
252  Real radii[3] = {elem->point(0).norm(), elem->point(1).norm(), elem->point(2).norm()};
253 
254  // Compute absolute differences between radii so we can determine which two are on the same
255  // circular arc.
256  Real dr[3] = {std::abs(radii[0] - radii[1]),
257  std::abs(radii[1] - radii[2]),
258  std::abs(radii[2] - radii[0])};
259 
260  // Compute index of minimum dr.
261  auto index = std::distance(std::begin(dr), std::min_element(std::begin(dr), std::end(dr)));
262 
263  // Make sure that the minimum found is also (almost) zero.
264  if (dr[index] > TOLERANCE)
265  mooseError("Error: element had no sides with nodes on same radius.");
266 
267  // Get list of all local node ids on this side. The first
268  // two entries in nos correspond to the vertices, the last
269  // entry corresponds to the mid-edge node.
270  nos = elem->nodes_on_side(index);
271 
272  // Compute the angles associated with nodes nos[0] and nos[1].
273  Real theta0 = std::atan2(elem->point(nos[0])(1), elem->point(nos[0])(0)),
274  theta1 = std::atan2(elem->point(nos[1])(1), elem->point(nos[1])(0));
275 
276  // atan2 returns values in the range (-pi, pi). If theta0
277  // and theta1 have the same sign, we can simply average them
278  // to get half of the acute angle between them. On the other
279  // hand, if theta0 and theta1 are of opposite sign _and_ both
280  // are larger than pi/2, we need to add 2*pi when averaging,
281  // otherwise we will get half of the _obtuse_ angle between
282  // them, and the point will flip to the other side of the
283  // circle (see below).
284  Real new_theta = 0.5 * (theta0 + theta1);
285 
286  // It should not be possible for both:
287  // 1.) |theta0| > pi/2, and
288  // 2.) |theta1| < pi/2
289  // as this would not be a well-formed element.
290  if ((theta0 * theta1 < 0) && (std::abs(theta0) > 0.5 * libMesh::pi) &&
291  (std::abs(theta1) > 0.5 * libMesh::pi))
292  new_theta = 0.5 * (theta0 + theta1 + 2 * libMesh::pi);
293 
294  // The new radius will be the radius of point nos[0] or nos[1] (they are the same!).
295  Real new_r = elem->point(nos[0]).norm();
296 
297  // Finally, move the point to its new location.
298  elem->point(nos[2]) = Point(new_r * std::cos(new_theta), new_r * std::sin(new_theta), 0.);
299  }
300  }
301 
302  return dynamic_pointer_cast<MeshBase>(mesh);
303 }
void addRequiredRangeCheckedParam(const std::string &name, const std::string &parsed_function, const std::string &doc_string)
These methods add an range checked parameters.
MetaPhysicL::DualNumber< T, D > abs(const MetaPhysicL::DualNumber< T, D > &in)
std::unique_ptr< MeshBase > generate() override
Generate / modify the mesh.
const Real _outer_radius
Radius of the outer circle. Logically, it&#39;s bigger that inner_radius.
const boundary_id_type _cylinder_bid
The boundary id to use for the cylinder.
Real _radial_bias
Factor to increase initial_delta_r for each ring.
const bool _use_tri6
Generate mesh of TRI6 elements instead of TRI3 elements.
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
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.
void mooseError(Args &&... args) const
Definition: MooseObject.h:147
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
unsigned int _num_rings
Number of rings.You can&#39;t specify both the number of rings and the radial bias if you want to match a...
SpiralAnnularMeshGenerator(const InputParameters &parameters)
registerMooseObject("MooseApp", SpiralAnnularMeshGenerator)
Real pow(Real x, int e)
Definition: MathUtils.C:211
Generates a spiral annular mesh given all the parameters.
InputParameters validParams< SpiralAnnularMeshGenerator >()
const Real _inner_radius
Radius of the inner circle.
PetscInt n
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...
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...
const unsigned int _nodes_per_ring
Number of nodes on each ring.
InputParameters validParams< MeshGenerator >()
Definition: MeshGenerator.C:16
const boundary_id_type _exterior_bid
MeshGenerators are objects that can modify or add to an existing mesh.
Definition: MeshGenerator.h:30