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 
19 {
22  "inner_radius", "inner_radius>0.", "The size of the inner circle.");
23  params.addRequiredRangeCheckedParam<Real>("outer_radius",
24  "outer_radius>0.",
25  "The size of the outer circle."
26  " Logically, it has to be greater than inner_radius");
27  params.addRequiredRangeCheckedParam<unsigned int>(
28  "nodes_per_ring", "nodes_per_ring>5", "Number of nodes on each ring.");
29  params.addParam<bool>(
30  "use_tri6", false, "Generate mesh of TRI6 elements instead of TRI3 elements.");
31  params.addRequiredRangeCheckedParam<unsigned int>(
32  "num_rings", "num_rings>1", "The number of rings.");
33  params.addParam<boundary_id_type>(
34  "cylinder_bid", 1, "The boundary id to use for the cylinder (inner circle)");
35  params.addParam<boundary_id_type>(
36  "exterior_bid", 2, "The boundary id to use for the exterior (outer circle)");
37  params.addParam<Real>("initial_delta_r",
38  "Width of the initial layer of elements around the cylinder."
39  "This number should be approximately"
40  " 2 * pi * inner_radius / nodes_per_ring to ensure that the"
41  " initial layer of elements is almost equilateral");
42  params.addClassDescription("Creates an annual mesh based on TRI3 elements"
43  " (it can also be TRI6 elements) on several rings.");
44 
45  return params;
46 }
47 
49  : MooseMesh(parameters),
50  _inner_radius(getParam<Real>("inner_radius")),
51  _outer_radius(getParam<Real>("outer_radius")),
52  _radial_bias(1.0),
53  _nodes_per_ring(getParam<unsigned int>("nodes_per_ring")),
54  _use_tri6(getParam<bool>("use_tri6")),
55  _num_rings(getParam<unsigned int>("num_rings")),
56  _cylinder_bid(getParam<boundary_id_type>("cylinder_bid")),
57  _exterior_bid(getParam<boundary_id_type>("exterior_bid")),
58  _initial_delta_r(2 * libMesh::pi * _inner_radius / _nodes_per_ring)
59 {
60  // catch likely user errors
62  mooseError("SpiralAnnularMesh: outer_radius must be greater than inner_radius");
63 }
64 
65 std::unique_ptr<MooseMesh>
67 {
68  return std::make_unique<SpiralAnnularMesh>(*this);
69 }
70 
71 void
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  {
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  BoundaryInfo & boundary_info = mesh.get_boundary_info();
125 
126  // Data structure that holds pointers to the Nodes of each ring.
127  std::vector<std::vector<Node *>> ring_nodes(_num_rings);
128 
129  // Initialize radius and delta_r variables.
131  Real delta_r = _initial_delta_r;
132 
133  // Node id counter.
134  unsigned int current_node_id = 0;
135 
136  for (std::size_t r = 0; r < _num_rings; ++r)
137  {
138  ring_nodes[r].resize(_nodes_per_ring);
139 
140  // Add nodes starting from either theta=0 or theta=pi/nodes_per_ring
141  Real theta = r % 2 == 0 ? 0 : (libMesh::pi / _nodes_per_ring);
142  for (std::size_t n = 0; n < _nodes_per_ring; ++n)
143  {
144  ring_nodes[r][n] = mesh.add_point(Point(radius * std::cos(theta), radius * std::sin(theta)),
145  current_node_id++);
146  // Update angle
147  theta += 2 * libMesh::pi / _nodes_per_ring;
148  }
149 
150  // Go to next ring
151  radius += delta_r;
152  delta_r *= _radial_bias;
153  }
154 
155  // Add elements
156  for (std::size_t r = 0; r < _num_rings - 1; ++r)
157  {
158  // even -> odd ring
159  if (r % 2 == 0)
160  {
161  // Inner ring (n, n*, n+1)
162  // Starred indices refer to nodes on the "outer" ring of this pair.
163  for (std::size_t n = 0; n < _nodes_per_ring; ++n)
164  {
165  // Wrap around
166  unsigned int np1 = (n == _nodes_per_ring - 1) ? 0 : n + 1;
167  Elem * elem = mesh.add_elem(new Tri3);
168  elem->set_node(0) = ring_nodes[r][n];
169  elem->set_node(1) = ring_nodes[r + 1][n];
170  elem->set_node(2) = ring_nodes[r][np1];
171 
172  // Add interior faces to 'cylinder' sideset if we are on ring 0.
173  if (r == 0)
174  boundary_info.add_side(elem->id(), /*side=*/2, _cylinder_bid);
175  }
176 
177  // Outer ring (n*, n+1*, n+1)
178  for (std::size_t n = 0; n < _nodes_per_ring; ++n)
179  {
180  // Wrap around
181  unsigned int np1 = (n == _nodes_per_ring - 1) ? 0 : n + 1;
182  Elem * elem = mesh.add_elem(new Tri3);
183  elem->set_node(0) = ring_nodes[r + 1][n];
184  elem->set_node(1) = ring_nodes[r + 1][np1];
185  elem->set_node(2) = ring_nodes[r][np1];
186 
187  // Add exterior faces to 'exterior' sideset if we're on the last ring.
188  // Note: this code appears in two places since we could end on either an even or odd ring.
189  if (r == _num_rings - 2)
190  boundary_info.add_side(elem->id(), /*side=*/0, _exterior_bid);
191  }
192  }
193  else
194  {
195  // odd -> even ring
196  // Inner ring (n, n+1*, n+1)
197  for (std::size_t n = 0; n < _nodes_per_ring; ++n)
198  {
199  // Wrap around
200  unsigned int np1 = (n == _nodes_per_ring - 1) ? 0 : n + 1;
201  Elem * elem = mesh.add_elem(new Tri3);
202  elem->set_node(0) = ring_nodes[r][n];
203  elem->set_node(1) = ring_nodes[r + 1][np1];
204  elem->set_node(2) = ring_nodes[r][np1];
205  }
206 
207  // Outer ring (n*, n+1*, n)
208  for (std::size_t n = 0; n < _nodes_per_ring; ++n)
209  {
210  // Wrap around
211  unsigned int np1 = (n == _nodes_per_ring - 1) ? 0 : n + 1;
212  Elem * elem = mesh.add_elem(new Tri3);
213  elem->set_node(0) = ring_nodes[r + 1][n];
214  elem->set_node(1) = ring_nodes[r + 1][np1];
215  elem->set_node(2) = ring_nodes[r][n];
216 
217  // Add exterior faces to 'exterior' sideset if we're on the last ring.
218  if (r == _num_rings - 2)
219  boundary_info.add_side(elem->id(), /*side=*/0, _exterior_bid);
220  }
221  }
222  }
223 
224  // Sanity check: make sure all elements have positive area. Note: we
225  // can't use elem->volume() for this, as that always returns a
226  // positive area regardless of the node ordering.
227  // We compute (p1-p0) \cross (p2-p0) and check that the z-component is positive.
228  for (const auto & elem : mesh.element_ptr_range())
229  {
230  Point cp = (elem->point(1) - elem->point(0)).cross(elem->point(2) - elem->point(0));
231  if (cp(2) < 0.)
232  mooseError("Invalid elem found with negative area");
233  }
234 
235  // Create sideset names.
236  boundary_info.sideset_name(_cylinder_bid) = "cylinder";
237  boundary_info.sideset_name(_exterior_bid) = "exterior";
238 
239  // Find neighbors, etc.
240  mesh.prepare_for_use();
241 
242  if (_use_tri6)
243  {
244  mesh.all_second_order(/*full_ordered=*/true);
245  std::vector<unsigned int> nos;
246 
247  // Loop over the elements, moving mid-edge nodes onto the
248  // nearest radius as applicable. For each element, exactly one
249  // edge should lie on the same radius, so we move only that
250  // mid-edge node.
251  for (const auto & elem : mesh.element_ptr_range())
252  {
253  // Make sure we are dealing only with triangles
254  libmesh_assert(elem->n_vertices() == 3);
255 
256  // Compute vertex radii
257  Real radii[3] = {elem->point(0).norm(), elem->point(1).norm(), elem->point(2).norm()};
258 
259  // Compute absolute differences between radii so we can determine which two are on the same
260  // circular arc.
261  Real dr[3] = {std::abs(radii[0] - radii[1]),
262  std::abs(radii[1] - radii[2]),
263  std::abs(radii[2] - radii[0])};
264 
265  // Compute index of minimum dr.
266  auto index = std::distance(std::begin(dr), std::min_element(std::begin(dr), std::end(dr)));
267 
268  // Make sure that the minimum found is also (almost) zero.
269  if (dr[index] > TOLERANCE)
270  mooseError("Error: element had no sides with nodes on same radius.");
271 
272  // Get list of all local node ids on this side. The first
273  // two entries in nos correspond to the vertices, the last
274  // entry corresponds to the mid-edge node.
275  nos = elem->nodes_on_side(index);
276 
277  // Compute the angles associated with nodes nos[0] and nos[1].
278  Real theta0 = std::atan2(elem->point(nos[0])(1), elem->point(nos[0])(0)),
279  theta1 = std::atan2(elem->point(nos[1])(1), elem->point(nos[1])(0));
280 
281  // atan2 returns values in the range (-pi, pi). If theta0
282  // and theta1 have the same sign, we can simply average them
283  // to get half of the acute angle between them. On the other
284  // hand, if theta0 and theta1 are of opposite sign _and_ both
285  // are larger than pi/2, we need to add 2*pi when averaging,
286  // otherwise we will get half of the _obtuse_ angle between
287  // them, and the point will flip to the other side of the
288  // circle (see below).
289  Real new_theta = 0.5 * (theta0 + theta1);
290 
291  // It should not be possible for both:
292  // 1.) |theta0| > pi/2, and
293  // 2.) |theta1| < pi/2
294  // as this would not be a well-formed element.
295  if ((theta0 * theta1 < 0) && (std::abs(theta0) > 0.5 * libMesh::pi) &&
296  (std::abs(theta1) > 0.5 * libMesh::pi))
297  new_theta = 0.5 * (theta0 + theta1 + 2 * libMesh::pi);
298 
299  // The new radius will be the radius of point nos[0] or nos[1] (they are the same!).
300  Real new_r = elem->point(nos[0]).norm();
301 
302  // Finally, move the point to its new location.
303  elem->point(nos[2]) = Point(new_r * std::cos(new_theta), new_r * std::sin(new_theta), 0.);
304  }
305  }
306 }
static InputParameters validParams()
Typical "Moose-style" constructor and copy constructor.
Definition: MooseMesh.C:78
void addRequiredRangeCheckedParam(const std::string &name, const std::string &parsed_function, const std::string &doc_string)
These methods add an range checked parameters.
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template * sin(_arg) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(tan
registerMooseObject("MooseApp", SpiralAnnularMesh)
const Real radius
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.
MeshBase & mesh
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
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
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...
int8_t boundary_id_type
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template cos(_arg) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(cos
MeshBase & getMesh()
Accessor for the underlying libMesh Mesh object.
Definition: MooseMesh.C:3198
SpiralAnnularMesh(const InputParameters &parameters)
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
Definition: MooseMesh.h:88
libmesh_assert(ctx)
const Real _initial_delta_r
static InputParameters validParams()
virtual void buildMesh() override
Must be overridden by child classes.
const boundary_id_type _cylinder_bid
The boundary id to use for the cylinder.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void mooseError(Args &&... args) const
Emits an error prefixed with object name and type.
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:2849
const unsigned int _nodes_per_ring
Number of nodes on each ring.
MooseUnits pow(const MooseUnits &, int)
Definition: Units.C:537
void ErrorVector unsigned int
const Real pi