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