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