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